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I. Introduction 


The purpose of this documentation is to familiarize the user with the theory and numerical 
treatment used in the 2-level GLA global fourth-order shallow water model. This model was de- 
signed to emulate the horizontal finite-differences used by the GLA Fourth-Order General Circula- 
tion Model (see Kalnay et al., 1983) in addition to its grid structure, form of high-latitude and 
global filtering, and time-integration schemes. Although shallow water models have been previously 
used to compare with GCMs, they have the disadvantage that, being barotropic, they cannot 
reproduce one of the most fundamental processes in the atmosphere. For this reason we have decid- 
ed to develop a 2-layer shallow water model, which can also be used as a conventional 1 -layer 
model. 

It has been found that this model provides a very useful tool in assessing the impact of 
numerical techniques under consideration for use in the GLA GCM. Experiments on the high- 
latitude and global filter (Takacs and Balgovind, 1983a, b), linear and non-linear normal mode in- 
itialization (Navon et al., 1985, Takacs, et al. 1985), aerosol tracer advection, and multi-variate data 
assimilation have been performed using this model with satisfying results. 

Section II develops the equations used for the two-layer shallow water system. Section III 
describes the grid structure and the horizontal finite-difference equations used by the model, while 
Section IV considers the addition of friction and external forcing to the model. The time-integration 
schemes used by the model are presented in Section V. In Section VI an examination is made of the 
types of filters needed by the model, i.e., high-latitude filtering to control linear instability due to 
fast-moving gravity waves near the poles, and global filtering to control non-linear instability due to 
aliasing. Section VII presents the normal mode initialization process and its implementation into the 
model, while Section VIII discusses the addition of a passive tracer. In Section IX a user’s guide is 
provided instructing the user on how to create intitial conditions, execute the shallow water model. 
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and post-process the data history. Section X contains the list of references used, while Section XI 
contains sample outputs. 

II. Model Equations 
a. Away from pole 

The governing equations for the shallow water model are based on a two-layer inviscid fluid 
system using the quasi-static approximation. Referring to Fig. 2.1, the momentum and hydrostatic 
equations for each layer are given by 



Fig. 2.1 


, d a , v „ 
where dTm + x > , ' v - 

Integrating (2.2) from any given height z within the layer to the top of the model, we have 
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Layer 1: 


l ^ dz=-plg ( hT . -z )> z - hx ’ ( 2 - 3 ) 

Layer 2: 

f 2 ^dz + | ^dz= -p 2 g(h T2 -z)-p 1 g(h Tl -h T2 ), z<h T2 (2.4) 

From (2.3) and (2.4) we see that 

P 2 (z) = pig(h Tl - z), z > h T2 (2.5) 

and 

P 2 (z) = p 2 g(h T2 - z) + Pig(h Tl - h T2 ), z < h T2 (2.6) 

Here we used Pi(h Tl ) = 0 and Pi(h X2 ) = P 2 (h T2 ). Taking the horizontal gradient of (2.5) and (2.6) 
yields 

-J-VP,=gVh Tl (2.7) 

Pi 

and 

-J- VP 2 =gVh T2 + £ gV(h Tl - h T2 ). (2.8) 

P2 p2 

Using (2.7) and (2.8), the momentum equations become 

^ Vi=-f kx Vi-gVhx, (2.9) 

and 

— \y 2 = - f k x \y 2 - gVhx 2 - gV(h Tl - h T2 ). (2.10) 

or, in general for layer k, 

^ \y k = — f k x \y k — gV[a k h T , + (1 — a k )h T2 ] (2.11) 
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where a x = 1, a 2 = pi/p 2 , and 


D d 

Dt k at +Nyk ' v - 


The continuity equation for layer k is given by 


„ aw k 

v-^ + ir =0 - 


Integrating (2.12) from h Tk+ , to h Tk yields 


dh k 

^=-v-(h k v k ) 


In summary, the momentum and continuity equations for layer k are given by 
Vk= -(Vk- V)\y k -f kx \y k -gV[a k h Tl + (l-a k )h T2 ] 




The flux form of the shallow water equations may be obtained by multiplying the momentum 
equation by h k and adding it to the continuity equation multiplied by V k . Doing so we find that the 
shallow water equations in component form on a sphere are given by 


d(hu) 

at 


^(hu)u + ^(hv cos 4>)u + |f + u ^j hv 


d(hv) 

at 


l _a 

a cos 4> a> 


a cos <j> al tahTl + V a ) hr J 
TcoT$ ^ hu ^ v + 4> ^ hv cos ^ 

_ ?4 cahT1 + (1-a)hTJ 


^ | u tan <|) 


ah _ 1 a(hu) a(hv cos cj>) 

at a cos <f> [ a\ acj> 
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(2.18) 



In equations (2. 16)— (2. 18), we have dropped the subscript k for simplicity, 
b. At the pole 

The equations used at the pole are obtained, following Kalnay et al. (1983), by first transform- 
ing the equations using a stereographic projection, and then area-averaging the equations over a 
polar cap. 

The stereographic projection of u and v result in a unique U p and V p , independent of longitude 


X. These are given by 

U pm = -u sin (X) - (-l) m v cos (X) (2.19) 

V p m = (-l) m u cos (X) - v sin (X) (2.20) 

where m = l for the south pole, and m=2 for the north pole. Multiplying (2.16) and (2.17) by the 
appropriate sin (X) or cos (X) and adding, we obtain 

l hp " =_ r^[l; (hu)+ ^ (hvcos<t>) ] < 2 - 21 ) 

Ft (hUp " ) = [s (hu)Up " + <hv cos <WVp "] + fhVp ” 

+ t[^S + <“ 1 >" cos ^] < ob " + <! - “) w < 2 - 22 ) 

Ft (hVpJ “ (hu) Vp " + 5? (hv cos ' t ' )Vp ”] " fhUp " 

- f [ ( - 1)m i - sin " 35] <“ h * + < 2 - 23 > 


Equations (2.21)-(2.23) are then area-averaged over a small polar cap. The results of this integration 
will be shown in Section III. 
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III. Horizontal finite differences 


a. Away from pole 

Following Kalnay et al. (1983), the GLA shallow water model uses a non-staggered A-grid 
(Arakawa, 1972), i.e. the height and wind components are defined at every grid-point: 



h,\y 

h,\y 

h,Ny 

1 

> . 
A<j> - 

h,\y 

h,\y 

h,\y 

1 

•” 

h,\y 

• ( 

AX 

h,\y 

h,\y 

X | 


> 1 

r ” 


i — 1 i i + 1 

Fig. 3.1 


Defining IM as the number of grid-points in the zonal direction and JM as the number of 
grid-points in the meridional direction, then the grid increments and cj>j are given by 


\ ( = -17 + (i - 1)A\, AX = 2 tt/IM 

(3.1) 

4>j = — tt/ 2 + (j — 1)A4>, Ac}) = tt/(JM — 1) 

(3.2) 

As an aid in notation, let us define the following: 



(3.3) 

8 (U 2 X) _ qf^x - qf-i 

M q )j 2 Ax 

(3.4) 

— x _ Qj + 1 Qj 

Hj + V2 ~ 2 

(3.5) 

„2x _ q j+2 + q j 

q j+' 2 

(3.6) 
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Then * 7 -^ may be approximated to fourth-order accuracy by 
ox 


§) = 5 Mq')/ - 1 + o(ax 4 ). 


When 5 operates on a single quantity as in the above example, (3.7) simplifies to 


dq\ 4/qj+i — qj-A 1 /qj+z — qj-a 


dxj . 3\ 2Ax 


3 \ 4Ax 


+ o(Ax 4 ). 


Before writing the finite-difference version of the momentum and continuity equations, let us 


define the mass variables U and V as 


and the generalized height h T as 


U = hu 


V = hv cos 4> 


h T = ah-r, + (1 — a)h Tr 


Then, using the notation of (3.3)— (3.6) together with the definitions (3.9)— (3. 1 1), the model equa- 
tions (2.16)-(2.18) may be written in finite-difference form as 


b = a -ib [1 8 > <u ‘ - 5 wu “ + ! 8 ‘ (v * - 1 8 *< v! *' ,! *> 


+ ( f + ^ H, - ^ [I Wft) - 1 W») 

\ ' i , j J L J 

Tt (hv)i > “ a-^b [f M0V) - 1 W02Xv!>) + 5 8 ‘ (W) - 5 


, . u tan cj>\ ,, . 

f + — ^ (hu)i.j' 


’ 3 n 52<i>(hT > ) 


l ^ = rsb t - 1 8 “ (ua) + 1 8 *< v *> - 1 • 
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It can be shown that (3.12)-(3.I4) conserves total energy in finite-difference form. Noting that 


the time rate-of-change of kinetic and potential energy is given by 


lhi(u= + v* + gh, + gh.) = u%> + + v>)f + gh T fJ (3. 


. 15 ) 


where h s is the height of the surface topography, then after some manipulation we find 


-2 

at .. 

>.j 


l 


= o 


hi, j 2 (u 2 + v 2 + gh T + gh s )i,j a 2 cos cfo AX A<j> 
which is the finite-difference analog of 

,2-n ,t!2 * 

— h^(u 2 + v 2 + gh T + gh s ) a 2 cos <J> dX d<> = 0. 


( 3 . 16 ) 


(3.17) 


b. At the pole 

As previously mentioned in Section 2, the equations used at the poles are obtained by first 
transforming the governing model equations using a stereographic projection, and then integrating or 
area-averaging these equations over a small polar cap. 


Let us consider integrating some function F(\, <f>) over a north polar cap: 

r 2i t -tt/2 


Fnp = [ [ F(\, (j)) a 2 cos c}) d\ dc)>. 

J 0 A<|> 


Similarly, for the south pole we have 


J p 27 T /• — Tt /2 + A<}) 

F (X, 4>) a 2 cos (j) dX d(j>. 

0 


( 3 . 18 ) 


( 3 . 19 ) 


Letting 4>' = -<j), (3.19) becomes 


n -rr/2 

F(X, -(j) 1 ) a 2 cos (f) 1 dX dcf) 1 . 

t/2— A<|> 


(3.20) 


Since 4> and c}) 1 are dummy variables in (3.19) and (3.20), we see that the only difference bet- 
ween integrating F (X, <t>) over the north pole and over the south pole is the substitution of — 4> 
for in the function F (X, (J>). 
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We shall now consider integrating the transformed equations (2.21)-(2.23). Let us define as 
the area of the polar cap, given by 


J -2tt /-tt/2 

I a 2 cos 4> d\ d<J> = 2-iTa 2 (1 — cos Ac|>). 

0 •'ir/2-A* 


(3.21) 


Then, integrating the continuity equation, we have 


ah P = _J_ f 
at A c J 0 


r2-n 

rir/2 

I • 

Kf2-^ 


+ (-l) m ^ ( hv cos 4>) a d\ d<t> 


(3.22) 


where m = 1 for the south pole and m = 2 for the north pole. The first term on the RHS of (3.22) 
vanishes since it is being integrated around a complete latitude circle. The second term, however, 
becomes 

-f f (-l) m tt (hv cos <j>) a dX d<f> = (-l) m a f (hv cos <j>) / dX . (3.23) 

Jo JH/2-A4. a<P J 0 /-IT/2-A4. 

Equation (3.23) is approximated in finite-difference form as 

(-l) m a f 2lT (hv cos c|>)/dX - (-l) m 27Ta T S 1 !? — 2 (hv),.! po ,e (3.24) 

Jo 'lr/2-A<t> 1M i = l 

Substituting (3.24) into (3.22) we have 

= <-‘> m a IM(i n - A cts -A40 I (hV) ' ' - - (3 ' 25) 

In practice, in order to achieve fourth-order accuracy, we repeat the process over a polar cap of 
2 A <f>, yielding 


t_ r-n m sin 2A<j> y rhvi 
at P ( * a IM (1 - cos 2A4>) ^ ^ ^ 


pole ±2- 


Combining (3.25) and (3.26) we obtain 


A<fr 2A<t> 

ah 4 ah l ah 

at p 3 at p 3 at p 
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This entire process is also applied to the transformed momentum equations (2.22) and (2.23), 


yielding 


|(hUp) = H(hUrt i ‘-i|;(hUp) 


2A< t> 


3 at’ 


3 at' 


(3.28) 


and 


4 a 


m ia 


-<hv P )=^ ( hv P ) -~(hvp) 


-2A4> 


(3.29) 


where 


_nA<J> 


at 


(hUp) = C(nA4>) 


IM 


i— 1 


(hvUp m ), J pole ±n 

ghpole ^Tj,j pole ± n COS \j 


IM 


i=l 


fpole hpole Vp m (3.30) 


_nA<£> 


^(hVp) = C(nAcJ)) 


IM 


(-l) m 2 (hvVp m )i > j pole ± n 


and 


( S^pole ^Tj j pole ± n Sin \j 

i=l 

sin (nAcf)) 


fpole hpole U Pm (3.31) 


C(nA4>) = 


a IM (1 — cos (nA<j>)) 
h T = ah T , 4" (1 ot) h T2 . 


IV. Friction and forcing 

Frictional damping and zonal forcing can be added as source and sink terms to the momentum 
equations as input through the NAMELIST parameters. To examine this, let us consider the tenden- 
cy of the zonal wind component, neglecting all other processes, as 


3u _ (u ~ u) 

at t 


(4.1) 
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where r is a time constant, and n is the mean zonal wind as a function of latitude. Integrating (4.1) 
in time from t to t + At, we find 


u(t + At) = n(l - e _At/T ) + u(t)e _At/T . 

(4.2) 

Thus, at At = 0 we have 

u = u(t), 

(4.3) 

while for At— we obtain 

u = u. 

(4.4) 

Letting a = At h, we may expand the exponential as 

e“ = 1 — a + a 2 — a 3 + ••• 

(4.5) 

Using (4.5) in (4.4) we find 

u(t + At) = u(t) + (u — u(t)) (a — a 2 + a 3 — • • •) . 

(4.6) 

In practice, (4.6) is truncated after the quadratic term, giving 

u(t + At) = u(t) + a(l — a) (u — u(t)) . 

(4.7) 


The frictional damping and zonal forcing terms are added after all other hydrodynamic pro- 
cesses are completed. During the model integration, both the u and v wind components are updated, 
although the zonal mean v wind is assumed to be zero. The mean zonal u wind is computed from 
the initial condition. The dissipative time constant r is a NAMELIST parameter. 

V. Time integration schemes 

Following Kalnay et al. (1983), the model has the option of using either the Matsuno time 
differencing scheme, the leapfrog time scheme, or a combination of both. The choice of which 
scheme is used is determined by the SWE NAMELIST parameter SCHEME which is set to either 
‘MATS’ or ‘LEAP’. The NAMELIST parameter NLEAP is the number of leapfrog time-steps to 
perform before restarting with a Matsuno time-step. This is done in order to control the separation 
of solutions arising from the computational mode of the leapfrog scheme. 
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The Matsuno time scheme is implemented using a predictor-corrector sequence. Letting q n 
represent a typical variable which is to be updated to time level n+1, and D(q") represent the non- 
linear space differences evalutated at time n, we have 

q* = q" + At D(q n ) (Predictor step) 
q n+1 = q" + At D(q*) (Corrector step) 

For the leapfrog time differencing scheme, we begin the forecast with the Matsuno scheme thus 
producing the n+1 time level, and then proceed with the leapfrog scheme using 

q n+ i = q n l + 2 At D(q n ) (Leapfrog) 

VI. Filters 

a) High-latitude 

The process of filtering in global grid-point models is fairly common in meteorological studies. 
High-latitude filtering is done in order to avoid the use of prohibitively short time-steps required due 
to fast-moving inertia-gravity waves near the poles. The specific filtering response needed to ensure 
stability is defined by the space and time differencing schemes used, and by the type of grid struc- 
ture in the model (Arakawa and Lamb, 1977). 

Takacs and Balgovind (1983a) have examined in detail various filtering techniques for use in 
this model. This work and further studies have shown that, using realistic initial conditions and 
topography, filtering the time-tendencies of the prognostic fields yields the closest simulation to a 
control run using a small time-step and no high-latitude filter. Other filtering techniques examined 
were the prognostic field filter, the energy-conserving filter, and the P.G. filter (see Takacs and 
Balgovind, 1983a). 

In all techniques used, the filtered quantities were fast-Fourier transformed into their wave con- 
stituents whose amplitudes were then altered by a wave-dependent damping function. In order to ob- 
tain the functional form of the damping function for the time-tendency filter, we shall consider the 
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linearized zonal momentum and continuity equations with a mean zonal wind U and a mean depth 
H. Using a non-staggered A grid (Arakawa, 1972) and fourth-order spatial differences we have 

I u,., = L-^X S00..I - m.T («• 




a cos 4>j 8 ^ ij 

8 ( h )i i) ' 

a cos <pj v /,, 7 

(6.1) 

a cos 4>j 8(h)i J 

H \ Fj 

a cos 4>j 8 ^ ,- 7 

(6.2) 


5 / \ _ 4 / £i + ij Qi-i.A _ l/ Qi + 2 .i 2 . j 

Sfqlu 3 1 2 AX I 3 I 4AX 


and ( ) Fi denotes that the quantity is filtered with filter function Fj. 

Assuming solutions of the form 

(u, h)i,j = Re {u j; hj exp (T[i0 - vt])} 

where 0 = kAX and I = V-l, and that ( ) Fi implies multiplying the Fourier amplitudes by F j5 
then (6.1) and (6.2) become 


. _ F / U _co_ A g co £ 

VU ’ ' \a cos 4>j AX U) a cos 4>j AX ’ 


£ — -p ( U co r . H co A 
V 1 1 \a cos 4>j AX 1 a cos cj>j AX U) 


where co = 4/3 sin 0 - 1/6 sin 20. Simplifying (6.5) and (6.6) we obtain 

(v - Ucoj) Uj - gcojhj = 0 
— HcOj Uj + (v — Ucoj) hj = 0 

where 

coFj 

Wj _ a cos cj>j AX ' 

For non-trivial solutions, it must be true that 

v = coj (U + Cg) 
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where 

C g = ViH 

For the leapfrog or Matsuno time integration schemes, linear stability requires 


Using (6.10) in (6.9) we obtain 


vAt < 1 . 


( 6 . 10 ) 


, ^ aAX cos <|)j 
j (U + Cg)At co 


In practice, we let 


Fj - mini 1, 


aAX cos 4>j 
(U + Cg)At co 


( 6 . 11 ) 


( 6 . 12 ) 


with U = 50 m/sec and C g = VgH.The mean height H is calculated using the initial conditions for 
a given experiment. Fig. 6.1 shows the damping coefficients as a function of wavenumber for the 
five latitudes nearest the poles (i.e. <{> = 86, 82, 78, 74, and 70 degrees) which are filtered in the 4 
deg. lat. x 5 deg. long, resolution model. Here we are using H = 10,000 m and At = 450 sec. 


High-Latitude Filter Response 
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We see that the longest and the shortest scales are not affected by the filter. The filter is the 
strongest near the 4 grid-length wave, which is characteristic of the fourth-order differences on the 
A grid. 

In summary, the time-tendencies of the prognostic fields, 

Fourier transformed into their wave consituents whose amplitudes are then multiplied by the damping 
function F. given by (6.12). These adjusted fields are then back-transformed and used to update the 
original fields. 

b. Global 

Spatial difference schemes which are not enstrophy conserving nor implicitly damping require 
global filtering of short waves to eliminate the build-up of energy in the shortest wavelengths due to 
aliasing. Global filtering as applied in the GLA GCM (see Kalnay et al., 1983) is performed by us- 
ing a 16th-order Shapiro (1979) filter on the sea-level pressure, sea-level temperature, potential 
temperature and winds. In the GLA shallow water model, the 16th-order Shapiro filter is applied to 
the wind and height fields at some specified frequency (normally every two simulated hours). 

The form of the two-dimensional 16th-order Shapiro filter is simply two passes of the one- 
dimensional filter applied in succession. To obtain the high order of the filter, eight passes of the 
3-point second-order operator are used as given by 

q, j = [l-(F|) 8 ][l-(FD 8 ]q i . j 

where 

F x(qc,j) = (qi+i.j - 2qi,j + qi_ lfj )/4 
F 4>(qi.j) (qi.j +1 2 q (> j + qj,j_i)/4 

Here, q is any quantity being filtered. 

To understand the effect of using the one-dimensional Shapiro filter on a given quantity, let us 
assume that 


(6.13) 

(6.14) 

(6.15) 


i.e. 


3h 

at 


a(hu) a(hv)\ 

at ’ at j 


are 


qj = Re {qe ikjAX } 


(6.16) 



Operating on qj using (6.14) yields 

F 2 (qj) = (qj+i - 2qj + q,_t)/4 (6.17) 

(e i9 - 2 + e -i9 ) 

= 4 qj (6.18) 

where 0 = kA\ 
or 

F 2 (qj) = -(sin 2 |J qj (6.19) 

Thus, the 16th-order Shapiro filter has a frequency response given by 

qj = (l-sin 16 ®)qe iie , (6.20) 

where each wavenumber amplitude coefficient q is modified by an amount (1 - sin 16 0/2). Fig. 6.2 
shows the damping characteristics of the 16th-order Shapiro filter as a function of wavenumber, in 
addition to the 2nd-, 4th- and 8th-order filters for comparison. We clearly see that the 2 grid-interval 
wave is damped completely after one application. Scales longer than 4 grid-lengths are almost un- 
touched for the 16th-order filter. 


Shapiro Filter 
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Takacs and Balgovind (1983b) have shown that filtering of the winds in the meridional direc- 
tion creates spurious sources of vorticity and divergence near the poles which ultimately affect the 
solution at all latitudes. Therefore, the 16th-order Shapiro filter is applied in practice to the height 
field in both directions and to the wind field in the zonal direction only. 

VII. Normal modes 

A number of experiments incorporating the normal modes of the GLA 4th-order shallow water 
model have been performed (see Navon et al. (1985), Takacs et al. (1985) ). These works are based 
on the work of Bloom (1983) in which the normal modes of the GLA 4th-order General Circulation 
Model were computed. 

The normal modes of the shallow water model have been incorporated into the shallow water 
system in two ways, a) Initialization, and b) High-Latitude Filtering. 

a. Initialization 

Following Daley (1980), based on the work of Machenhauer (1977), we can write the time- 
change of the normal mode projection coefficients C as 

dC . ^ 

dt -lw C + r > (7.1) 

where w is the normal mode eigenfrequency. Here r represents then non-linear and forcing terms. 
Assuming r to be independent of time, integrating (7.1) yields. 

c (t) = ( c (0 + i) e “-i. (7.2) 

The first term in (7.2) represents the oscillation term which we want to suppress for normal mode 
initialization. Clearly, if 

< 7 - 3 > 
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model integration step, and dividing by the normal mode eigenfrequency. This procedure is usually 
repeated until a critical number of iterations (4-10) have been performed, or some convergence 
criterion has been reached. 


It has been experimentally found that a few particular modes and wavenumbers fail to converge 
when the Machenhauer technique is applied, particularly mode 43 for wavenumber 1. However, a 
more robust method of initialization, based on the work of Rasch (1984), is also available in the 
shallow water system. In this method, the assumption of the non-linear terms in (7.1) being constant 
in time is no longer required. Defining 

F(C)-^, (7.7) 

then F may be expanded in a Taylor series about t=0, 

F(C(t)) = F(C c ) + (C - C 0 ) ^/ o , (7.8) 

or 


C 0 = C- 


F(C)-F(C 0 ) 
dF / 



(7.9) 


18 



By definition, however, we want F(Co) = 0, i.e. the initial normal mode coefficients whose 


time tendency for the fast modes are zero. Thus, (7.9) may be expressed as 

(C n +i - C n )/At 


c„ = c n - 


where we have approximated F(C) by 


F(C) : 


(dF/dC)/ c 


, C n+1 C n 

At 


(7.10) 


(7.11) 


Notice that from (7.1), ^ = iw if r is independent of C(t). Thus, under this assumption (7.10) 
reduces to the Machenhauer initialization previously discussed. In practice, however, we let 


where 


dF / _ F„ — F 0 
dC/ 0 C n - C 0 


(7.12) 


F 


n 



and 


F 0 - 



(7.13) 


(7.14) 


This procedure is again repeated until a critical number of iterations or a convergence criterion has 
been reached. 


The non-linear normal mode initialization techniques are summarized in Fig. 7.1. Here a plot 
is made of the time history of the LOG of the balance parameter BAL. BAL is defined to be the 
square of the time-tendencies of the normal mode projection coefficients, summed over all gravity 
modes and all wavenumbers. We clearly see that both the Machenhauer and Rasch techniques tend 
toward a balanced state during the first 5 iterations of the NMI scheme. After 5 iterations, however, 
the Machenhauer technique tends to diverge while the Rasch technique continues to reduce the 
balance parameter. It has been found that the failure by the Machenhauer technique is, as previously 
mentioned, primarily due to mode 43 for wavenumber 1, which results in a large divergent compo- 
nent in the polar wind field. 
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Non-Linear Normal Mode Initialization 



Rasch —& 
Mach. — * 


Iteration 

Fig. 7.1 


b. Filter 

Following the work of Navon et al. (1985) and Takacs et al. (1985), the model normal modes 
have been applied toward the problem associated with high-latitude filtering to control linear 
stability. 

In an explicit leapfrog or Matsuno time-integration scheme the maximum allowable time-step is 
governed by the frequency of the fastest propogating mode. For a given gravity mode with frequen- 
cy cog, the maximum permissible time-step is given by 

At max = l/|wg| . (7.15) 
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Thus, by ordering the normal modes in decreasing eigenfrequency, we may easily find all modes for 
which 

|tog| > 1/At (7.16) 

holds true. Here, At is the time-step used in the model integration. These modes, which consist of a 
relatively small number of gravity modes, are denoted as ‘fast modes’, while the remaining gravity 
modes and all of the Rossby modes are denoted as ‘slow modes’. Thus in order to prevent linear in- 

n 

stability, at each time-step we simply forward transform the dependent variables (u,v,h) given to us 
by the model thereby producing a normal mode coefficient set. To obtain the projection onto the fast 
modes, we reconstruct the fields by back-transforming only those mode coefficients whose eigenfre- 
quency is greater than the maximum allowed. This produces Au, Av, and Ah which are then sub- 
tracted from the current fields to produce ‘filtered’ new values. 

VIII. Tracer advection 

Through the use of the logical NAMELIST parameter LTRACE, it is possible to perform ex- 
periments to observe the time evolution of passive tracer advection. The numerical advection scheme 
used in these studies is presented in Takacs (1985), and has the property of adverting localized 
disturbances (or profiles) with minimized dissipation and dispersion errors. The actual version of the 
scheme implimented here is a slight variation from Takacs (1985) in the sense that both zonal and 
meridional advections are calculated simultaneously. 

The tracer advection scheme is derived from the conservation of the tracer having no source or 
sink terms. Letting q be defined as the mixing ratio of the tracer quantity, then 

4 / pqdV = f p^d V = f pSdV - f F • nds, (8.1) 

UtJ VoI(t) •'vol(t) at J VoJ(t) •'s 

where S is a source or sink of q per unit time, and F is the flux of q across the surface of the 
volume element. Equation (8.1) may be written as 

p^ = pS-V-F (8.2) 
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and for no sources or sinks and no net flux across the boundary, we have 


dq dq 


5q. 


dr = 7F +v 'Vq + w^=o 

For our shallow water system, we assume q = q(x,y,t), thus (8.3) becomes 

fjr + V-(\yq) = qV- \y. 

Using the shallow water continuity equation given by 

3h 


we see that 


— + \y -Vh + hV- \y =0, 


V- V =-I(f + V-V h ). 


Using (8.6) in (8.4) and re-arranging, we have 

a(hq) 


at 


+ V-(h\yq) = 0. 


Equation (8.7) may be written in component form on the sphere, given by 

a(hq) , 1 |"a(hu)q d 1 

— ^7 1 T rr f vr (hvCOS<£)q = 0. 

at a cos (j) L a\ d(j> ' M J 

Defining 

U = hu 
V = h v cos W>) 

8 X = a cos (^) 8\ 

8y - a8d> 

equation (8.8) may be written as 


a(hq) [ a(Uq) 1 1 a(v q ) 

at ax cosd> ay 


(8.3) 

(8.4) 

(8.5) 

( 8 . 6 ) 

(8.7) 

(8.8) 


(8.9) 
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Following Takacs (1985), equation (8.9) is approximated in a predictor-corrector sequence as 


follows: 


Defining 

^x O) = Fi + l/2,j — Fj— 1/2, j 

(G) == Gi,j + x/2 — Gj,j_i/2 

Ui+i/ 2 , j = (Ui+ij 4- Ui,j)At/Ax 
Vi, j+ i /2 = (Vi, j+1 + Vi.j) At/ Ay 

Then: 

PREDICTOR: 

(hqK, = (hq)r.,-[MF)+^S,(G)] 


Fi+l/2,j ^i + 1/2, j Qi, j + u i+i/2 f j ^r+ij 

Gi, j + 1/2 = V i, jH-1/2 Qi! j + Vi, J + l/2q"j + l 




i -H 1/2, j 


j H- 1/2 


CORRECTOR: 


i n+1 = n n — - 

iui qij 2[ 


s x (P) + 


cos^ 


,( R )] + 


M“Q) + 


— 8 y 
COS0J y 


OS) 


P i+ i/2,j = K + m,i (ql+i.j + q"j) + Ui+i/2.j (q*j + qP+i,j) 

Ri. j+1/2 = Vitj+1/2 (q*j+i + q"j) + Vi7 j+1/2 (q*j + q" j+1 ) 

Q|+W.J = Ui+i/2.j (qi+i.j - q"j - q*j + qP-i.i) + uf+i/2.] (q"+2,j ~ q*+i,j - qr + i,j + q*j) 

Si, j+1/2 = Vjj+,/ 2 (qitj+i - q"j - q*j + q" j-i) + v , ~ i+m (qr, j+2 - qfo*, - qS* >j+1 + qf j) 


Oti+i/2,j 


(i±M) 

' 6 'j+l/2,j 


Pi, j+1/2 — 
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IX. User’s guide 

The GLA Fourth-Order Shallow Water Model is run through use of a series of CMS EXEC’s 
kept on the F400M disk 200 machine. In order for the user to run the shallow water model, the user 
must be linked to a set of particular disks which contain model code, text and execs. Thus, the user 
must execute the following: 

LINKTO F400M 200 591 M 
LINKTO F400M 192 592 C 
LINKTO F400M 207 593 F 
590 3 

As previously mentioned, the F400M 200 disk contains the model code and CMS EXEC’s which 
run the model. The F400M 192 and 207 disks are library disks for the GLA modelling group which 
contain many useful execs and processing routines which are called upon from the shallow water 
utilities. The command 590 3 invokes 3 cylinders of TDISK at address 590 with filemode I. This is 
done so that all processing and output can be temporarily stored on TDISK and thereby not interfere 
with the user’s own A-disk. Should the user find that more than 3 cylinders are needed, this number 
can be increased. 

The shallow water exec system contains the following execs: SWEIC, SWERUN, SWEPOST, 
SWEPLOT, SWEBATCH, SWEMODE, SWEDIFF, SWEMVS. To obtain information on any of 
the execs, the user need only type the exec name followed by a question mark, eg. SWERUN ?A 
brief description of each exec and the namelist data set associated with each program is now 
presented. 
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SWEIC 


This program creates initial conditions for the GLA Fourth-Order 2-Level Shallow Water 
Model. The number of levels and the type of Initial Conditions desired are NAMELIST parameters. 
Output from this EXEC is the initial condition file SWEIC OUTPUT I. 

Hit return to continue, or X to quit. 

* ********************************************************************** 
** M/A-COM SIGMA DATA INC. NASA - GSFC * 
* ********************************************************************** 
* 

* THIS PROGRAM CREATES INITIAL CONDITIONS 

* FOR THE SHALLOW WATER MODEL 

* 

* ********************************************************************** 
* 

* INPUT : 

* 

* KMAX : Number of levels (1 OR 2) in model 

* SETUP : Choice of Initial Conditions 

* TPG : FLAG (I OR 0) for realistic topography 

* R : Wavenumber to be used for ROSSBY-HAURWI TZ wave 

* 

* LEVEL : Manditory pressure level(s) to be used 

* for realistic GCM Initial Conditions: 

* 


* 

Level = 

1 

50 

mb 

Level = 7 

300 

mb 

* 

Level = 

2 

70 

mb 

Level = 8 

400 

mb 

* 

Level = 

3 

100 

mb 

Level = 9 

500 

mb 

* 

Level = 

4 

150 

mb 

Level = 10 

700 

mb 

* 

Level = 

5 

200 

mb 

Level = 11 

850 

mb 

* 

Level = 

6 

250 

mb 

Level = 12 

1000 

mb 


* 


* NOTE: SETUP = 'RSBY' for ROSSBY-HAURWI TZ wave initial conditions 

* = 'GMCD' for UNINITIALIZED GCM DATA 

* = * WINT * for JAN 9, 1979 UNINITIALIZED DATA 

* = ' WNMI * for JAN 9, 1979 INITIALIZED DATA 

* = ' SUMM ' for JUNE 11, 1979 UNINITIALIZED DATA 

* = ' SNMI * for JUNE 11, 1979 INITIALIZED DATA 

* 

&LEVELS 
KMAX= 1 
&END 

&IC0N 

SETUP=’WINT’ , LEVEL=7 , 9 
&END 

&WAVEN0 
R=3, TPG=0 
&END 
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SWERUN 


This program runs the GLA Fourth-Order 2-Level Shallow Water Model. Experiment 
parameters such as the timestep, length of integration, frequency of output, etc. are NAMELIST 
parameters. The model reads initial conditions from the file SWEIC OUTPUT I. Output from the 
model is stored in the file VARS OUTPUT I. 


Hit return to continue, or X to quit. 


* ********************************* ******** ********** ************* ****** 

* * M/A-COM SIGMA DATA INC. NASA - GSFC * 

* 

* THIS PROGRAM RUNS THE GLA FOURTH-ORDER 

* 2-LEVEL SHALLOW WATER MODEL 

* 

* ********************************************************************** 
* 


* INPUT : 

* 

* SCHEME: 

* 

* NLEAP : 

* 

* TSTART : 

* 

* TMAX : 

* DT : 

* TPRT : 

* KAPPA : 

* KMAX : 

* 

* GFILT : 

* HFILT : 

* LNMF : 

* LNMI : 

* LRASCH: 

* NMIMAX : 

* TSHAP : 

* 

* LTRACE: 

* LPRIME : 

* LFORCE: 

* ALPHA : 

* 

* LFRICT: 

* TAO : 

* 

* XLAB : 


Time-scheme used in model, 'MATS' for MATSUNO 

'LEAP' for LEAPFROG 
Number of consecutive LEAPFROGS to do 
before restarting with a MATSUNO step 
Hour from INITIAL CONDITION 
at which to start integration 
Hour at which to stop integration 
Time-step in seconds 

Frequency (IN HOURS) at which to print output 
Ratio of densities between UPPER and LOWER level 
Number of levels used in model (1 or 2) 

Global Shapiro filter flag (T or F) 

High-latitude filter flag (T or F) 

High-latitude normal mode filter flag (T or F) 

Non-linear normal mode initialization flag (T or F) 

True for RASCH'S NMI iteration. False for MACHENHAUR 
Number of iterations for normal mode initialization 
Frequency (IN HOURS) at which to apply global shapiro filter 

Logical for advection of tracer 

Logical for experiments to remove initial time tendencies 

Logical for experiments to call FORCING terms 

Parameter to enhance or reduce FORCING ( alpha X FORCING ) 

Logical for frictional and zonal forcing 

Dissipation time scale (DAYS) for friction and zonal forcing 
Label associated with current experiment 
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&TSCHEM 

SCHEME = 'MATS’ , NLEAP=20 
&END 

&INPUT 

TSTART = 000.000, TMAX=000.50, DT=0450.0, TPRT=000.60 
&END 

&LEVELS 

KAPPA- 1.00, KMAX= 1 
&END 

&FILT 

GFILT= . T. , TSHAP=2 . , 

HFILT= . T . , LNMF =.F., 

LNMI = . T. , LRASCH=.T. , NMIMAX=10, 

&END 

&STRESS 

LFRI CT= . F . , TAO=7 . , 

&END 

&FORCE 

LPRIME= . F . , LFORCE= . F . , ALPHA = 50.0, 

LTRACE= . F . , 

&END 

&LABEL 

XLAB = ’RASCH TE',’CHNIQUE ’,’WITH NEW’,’ BALANCE’, 

’ PARAMET ' , ' ER 

&END 
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SWEPOST 


This program post-processes data from an SWERUN or SWEBATCH run. During the execu- 
tion of this EXEC, the input data is plotted and/or displayed in various user-defined formats. 
Diagnostic quantities such as vorticity, divergence, streamfunction, as well as zonal and eddy kinetic 
energies, are also available. All options are NAMELIST parameters. The input data set defaults to 
VARS OUTPUT I, but may be any stored output data set of the VARS OUTPUT format. 


Hit return to continue, or X to quit. 


* ************ ******************* ************** ************************* 
** M/A-COM SIGMA DATA INC. NASA - GSFC * 

* ********************************************************************** 
* 

* THIS PROGRAM POST-PROCESSES DATA FROM THE 

* 2-LEVEL SHALLOW WATER MODEL 

* 

* ************ ************************************************* ********* 
* 

* INPUT : 

* 

* ALL: Logical to read ALL data on history file 

* TIMIN: Beginning time for post -proce ss i ng if ALL = False 

* TIMOUT: Ending time for post-processing if ALL = False 

* 

* KU: FLAG for u wind component 

* KV : FLAG for v wind component 

* KH : FLAG for height field 

* KQ : FLAG for tracer 

* KVOR: FLAG for vorticity 

* KDIV: FLAG for mass divergence 

* KSTR: FLAG for streamfunction 

* KCHI : FLAG for velocity potential 

* KMDIV: FLAG for mass divergence 

* KPVOR: FLAG for potential vorticity 

* 

* GRF : FLAG for GRAFICS ROUTINE on CMS 

* DSP: FLAG for NUMERICAL DISPLAY of array 

* PHY:. FLAG for ENERGY and MOMENTUM analysis 

* ENS: FLAG for TIME SERIES of TOTAL POTENTIAL ENSTROPHY 

* and TOTAL MASS 

* 

* KFFT: Fourier decomposition of specified fields 

* 

* KPRS: FLAG to plot history record of sea-level pressure (heights) 

* NGRID: The total number of points used for KPRS 

* : IGRID: The i-index of points used for KPRS 

* JGRID: The j-index of points used for KPRS 

* SPLINE: FLAG for using SPLINE for KPRS plot 

* 
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* 

* SHAP: FLAG to Shapiro filter data before plotting 

* NSHAP: 1/2 the Order of the Shapiro filter (2,4,6, or 8) 

* DSHAP: The number of dimensions to filter, 1: zonal only, 2: both 

* 

* NUM: Offset in GRAFICS ROUTINE for UPPER and LOWER LIMITS (NUM=0,8) 

* NX: Number of interpolated points in x-direction for GRAFICS 

* NY: Number of interpolated points in y-direction for GRAFICS 

* Note: Usually NX=3, NY=2 for COARSE GRID 

* NX=2, NY=2 for FINE GRID 

* 

* NSCALE : FREQUENCY of scaling (NSCALE TIMES) 

* 

&PARMS 


ALL 

= F, 

TIMIN 

= 000., 

TIMOUT 

= 012. , 

GRF 

= F, 

DSP 

= F, 

PHY 

= F, 

ENS 

= F, 

KFFT 

= T, 

KQ 

= F, 

KU 

= F, 

KV 

= F, 

KH 

» T, 

KVOR 

= F, 

KDI V 

= F, 

KSTR 

= T, 

KCHI 

= T, 

KMDIV 

= F, 

KPVOR 

= F, 

SHAP = 

F, NSHAP = 2, DSHAP 

NX = 2 , 

NY=2, NUM=0, NSCALE=3 


&END 

&GRIDS 

KPRS = F, SPLINE = T, 

NGRID = 8, 

IGRID = 01, 10, 20, 30, 40, 50, 60, 72, 
JGRID = 01, 05, 10, 15, 20, 30, 40, 46, 
&END 
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SWEPLOT 


This program post-processes data from an SWERUN or SWEBATCH run. During the execu- 
tion of this EXEC, the input data is plotted using the WOLFPLOT plotting package. With this 
package, wind vectors and streamlines as well as diagnostic quantities such as vorticity, divergence, 
etc. are available. All options are NAMELIST parameters. This plots may be run inter-actively or 
under the BATCH facility. A NON-LABEL tape is required for plots using BATCH. 

Hit return to continue, or X to quit. 

Type B for BATCH job or 
C for CMS 

Type in the VMID and the NODES (max =3) of the MVS data set. 
eg.)W3LLT NODE1 NODE2 NODE3 

Enter the NONLABEL SCRATCH TAPE if 

Enter TIME in minutes. Default is 5 Minutes. 

Enter FILENAME for batch namelist. Default is SWEBATCH. 
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*********************************************************************** 

* M/A-COM SIGMA DATA INC. NASA - GSFC * 
*********************************************************************** 

* 

* THIS PROGRAM PLOTS DATA FROM THE 

* 2-LEVEL SHALLOW WATER MODEL 

* 

* ****************************************************** **************** 
* 

* INPUT : 


* 


* 

X 

LTPG : 

FLAG 

f or 

topography 

DTPG : 

Contour 

interval 

X 

* 

LU: 

FLAG 

for 

u wind component 

DU: 

Contour 

interval 

* 

LV: 

FLAG 

for 

v wind component 

DV: 

Contour 

interval 

* 

LH: 

FLAG 

for 

height field 

DH: 

Contour 

interval 

* 

LQ: 

FLAG 

for 

tracer 




* 

LVOR: 

FLAG 

for 

vort ic i ty 

DVOR: 

Contour 

interval 

* 

LDI V : 

FLAG 

for 

mass divergence 

DDI V : 

Contour 

interval 

* 

LSTR: 

FLAG 

for 

streamf unct i on 

DSTR : 

Contour 

interval 

* 

LCHI : 

FLAG 

f or 

velocity potential 

DCHI : 

Contour 

interval 

* 

LHDIV: 

FLAG 

f or 

mass divergence 

DHDIV 

: Contour 

interval 

* 

LPVOR : 

FLAG 

for 

potential vorticity 

DPVOR 

: Contour 

interval 

X 

* 

Note : 

Contour 

interval = 0 implies 

computer 

generated 

interval 


* 

* LSTRLN: FLAG for streamlines 

* LVEC: FLAG for wind vectors 

* VREF: Wind velocity whose magnitude is represented by 1 grid length 

* 

* ALL: FLAG to read and process ALL data in history file 

* TIMIN : Time (in hours) to begin processing data if ALL = .F. 

* TIMOUT: Time (in hours) to end processing data if ALL = .F. 

* NSTART: YYMMDD at which to begin for labeling purposes 

* 

* LAT1 : Southern latitudinal boundary 

* LAT2 : Northern latitudinal boundary 

* L0NG1: Western longitudinal boundary 

* L0NG2: Eastern longitudinal boundary 

* 

* LPRIME: FLAG to plot differences between history file 

* and the initial condition 

* 
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&INPUT 


LPRIME = F, 


LTPG 

= 

F, 

DTPG 

= 

200. , 

LVEC 

- 

T, 

VREF 


30.0, 

LSTRLN 

= 

F, 




LU 

= 

F, 

DU 

= 

0.0 

LV 

= 

F, 

DV 

= 

0.0 

LH 

= 

T, 

DH 

= 

060.0 

LQ 

= 

F, 

DQ 

= 

5.0 

LVOR 

= 

F, 

DVOR 

= 

0.0 

LDI V 

= 

F, 

DDI V 

= 

0.0 

LSTR 

= 

F, 

DSTR 

= 

0.0 

LCHI 

= 

F, 

DCHI 

= 

0.0 

LHDIV 

= 

F, 

DHDIV 

= 

0.0 

LPVOR 

= 

F, 

DPVOR 

= 

0.0 

ALL 

— 

T, 




TIMIN 

= 

000 




TIMOUT 

= 

000 

'• > 



NSTART 

= 

o. 





LAT1 = -090, LAT2 = 90, 

LONG 1 = -180, L0NG2 = 180, 

&END 


CreatingSWEPLOT BATCH I 

/JOB W3LLT SWEPLOT 
/QUERY PRINT=* 

/IDENT SWEPLT 
/ERROR PRINT=* 

/SET RATE = BASIC 
/SET DUMP = NO 
/SET NOTIFY = YES 
/SET TIME = 10 
/SET CARDS = 100000 
/SET LINES = 100000 
/SET SIZE = 2 
/SET TAPE = 1 

// SETUP TAPE = 123456 NL RI 181 
LINKTO F400M 200 254 M 

EXEC GETTAP 123456 NL RI 181 

EXEC SWEPLOTS 123456 PL0T1 TS0UT2 W3LLT TRACER DATA 
/* 


Do you wish to SUBMIT? 
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SWEBATCH 


This program runs the GLA Fourth-Order 2-Level Shallow Water Model under a BATCH en- 
vironment. Input NAMELISTS for SWEIC and SWERUN are created and stored as SWEIC 
BATCH# and SWE BATCH# files. The BATCH # is USER defined, with a default of 
BATCH. The data history file, VARS OUTPUT, will be sent back to the USER’S reader upon 
completion. 


Hit return to continue, or X to quit. 

Enter JOB #. Default is BLANK. 

Enter STORAGE SIZE #. Default is 2 Meg. 

Enter TIME in minutes. Default is 2 Minutes. 
Do you wish to compile a new version on BATCH? 

Creating SWEBATCH JOB I 


/JOB W3LLT SWEJOB 
/QUERY PRINT=* 

/I DENT SWE 
/ERROR PRI NT=* 

/SET RATE = PRIORITY 
/SET DUMP = NO 
/SET NOTIFY = YES 
/SET TIME = 2 
/SET CARDS = 100000 
/SET SIZE = 2 
CP LINK W3LLT 191 237 RR 
ACC 237 B/A 
EXEC SWEBM BATCH 
/* 


Do you wish toSUBMIT? 



SWEMODE 


This program analyzes data from an SWERUN or SWEBATCH run. During the execution of 
this EXEC, the input data is projected onto specific user-defined wavenumbers and normal modes, 
and then reconstructs the u, v and h field using only those wavenumbers and modes specified. The 
default input data set is VARS OUTPUT I, although the user can input any data set of the VARS 
OUTPUT format currently stored. The output data set is called MODE DATA I. Also, the normal 
mode amplitudes of the output data set as a function of wavenumber is produced in the file AMPL 
OUTPUT I. 


Hit return to continue, or X to quit. 


* * M/A-COM SIGMA DATA INC. NASA - GSFC * 

* ****************** Is******** ***************** ft************ It************ 

* * 

* This program PROJECTS data U,V and H onto the * 

* Shallow Water Model NORMAL MODES, and then * 

* reconstructs the data using ONLY user-defined * 

* WAVENUMBERS (00-36) and MODES (01-66,67,68). * 

* * 

* Note: For ALL wavenumbers and/or ALL modes, * 

* set logicals ALLWAV and/or ALLMOD to True. * 

* Otherwise, set to False and PICK individual * 

* wavenumbers and modes. * 

* * 
* ********************************************************************** 
* 

&INPUT 

ALLWAV = .F. , 

ALLMOD = .F., 

NWAVE = 1,2, 3, 4, 

NMODE = 45,46,47,52, 

HEADER = 'NO FILTE','R - TEND' , ' ENCY FIL' , ' TER 


&END 


Type in FILENAME FILETYPE FILEMODE of the data to be processed. 
Default is VARS OUTPUT I, or ... 

Type 7 if the data is on the Amdahl V7. 
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SWEDIFF 


This program differences two stored data sets, and produces the differenced u, v and h fields in 
the file VARS OUTPUT I. The input data sets may by stored data on CMS or MVS. In addition to 
the VARS OUTPUT file, the STATS between the two input fields (i.e., RMS, correlation coeffi- 
cient, etc.) are created in the file STATS OUTPUT I for the globe and various regional areas. One 
region is USER defined, given as a NAMELIST input data file called SWEDIFF DATA A. 


Hit return to continue, or X to quit. 


Type inFILENAME FILETYPE FILEMODE of the first data set to read, 
or ... 

Type7 

followed by theVMIDand theNODES(max=3) of theMVSdata set 
if the data is on the Amdahl V7. 
eg.) 7 W3ABC NODE 1 N0DE2 NODE3 


Type inFILENAME FILETYPE FILEMODE of the second data set to read, 
or ... 

Type7 

followed by theVMIDand theNODES(max=3) of theMVSdata set 
if the data is on the Amdahl V7. 
eg . ) 7 W3ABC NODE1 NODE2 NODE3 
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SWEMVS 


This program stores data from an SWERUN or SWEBATCH run onto Mass Storage. During 
the execution of this EXEC, the data is temporarily converted to a formatted LRECL 80 data set 
called VARS1 OUTPUT I, and is then sent to Mass Storage via MVS. The default name of the in- 
put data set is VARS OUTPUT I, although the user can input any data set of the VARS OUTPUT 
format currently stored. The name of the stored data set on MVS is of the form: 

VMID node 1 node2 node3 
eg) W3ABC SWE TEST DATA 

Hit return to continue, or X to quit. 
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XI. Appendix A 

a. Sample output 

The following is a brief sample of the type of initial conditions available, the realistic 
topography data set used, and a few sample runs. 
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Figure 2. Topography used in conjunction with realistic initial conditions 
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Fig. 7. Rossby-Haurwitz wave initial condition using 36x23 resolution 
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Fig. 8. Rossby-Haurwitz wave after 4 days, using dt=900 seconds 
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Fig, 9. Rossby-Haurwitz wave after 8 days 
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Fig. 10. Rossby-Haurwitz wave after 12 days 
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Fig. 12. Tracer experiment after 1 day 
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b. SWE source code 


The following is the FORTRAN source code for the 2-layer shallow water model. 


C ********************************************************************** 
C* M/A-COM SIGMA DATA INC. NASA - GSFC * 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 


********************************************************************** 


PROGRAM SWE (INPUT, OUTPUT, 

$ UNITS = INPUT, 

$ UNIT6 = OUTPUT, 

$ UNIT45 = VARS, 

$ UNIT46 = ICSET ) 

********************************************************* 


**** **** 

**** THIS PROGRAM SOLVES A TWO-LAYER SHALLOW WATER **** 
**** MODEL WITH TOPOGRAPHY ON A SPHERE. THE NUMBER **** 
**** OF LAYERS USED IS A NAMELIST PARAMETER. **** 

**** **** 

**** SPATIAL DERIVATIVES ARE APPROXIMATED BY **** 

**** ENERGY-CONSERVING FOURTH-ORDER DIFFERENCES. **** 

**** THE MATSUNO OR LEAPFROG TIME-SCHEME MAY BE **** 

**** USED AND IS A NAMELIST PARAMETER. **** 

**** **** 

**** THE POLES ARE TREATED AS SINGULAR POINTS AND **** 
**** A SPECIAL SOLUTION IS USED THERE. **** 

**** **** 

**** A SIXTEENTH-ORDER SHAPIRO FILTER IS PERIOD- **** 

**** ICALLY APPLIED TO THE HEIGHTS IN BOTH **** 

**** DIRECTIONS AND TO THE WINDS IN THE ZONAL **** 

**** DIRECTION TO PREVENT NON-LINEAR INSTABILITY. **** 

**** HIGH-LATITUDE FOURIER FILTERING IS DONE TO **** 

**** PREVENT LINEAR INSTABILITY. **** 

**** **** 


********************************************************* 


NAMELIST 

NAMELIST 

NAMELIST 

NAMELIST 

NAMELIST 

NAMELIST 


/TSCHEM/ SCHEME, 
/INPUT / TSTART, 
/LEVELS/ KMAX , 
/FILT / GFILT , 
/STRESS/ LFRICT, 
/FORCE / LPRIME, 


NLEAP 

TMAX, DT, TPRT 
KAPPA 

HFILT, LNMF, LNMI , LRASCH, 
TAO 

LFORCE, LTRACE, ALPHA 


NMIMAX, 


TSHAP 


LOGICAL LPRIME, LFORCE, LFRICT, LTRACE 
LOGICAL HFILT, GFILT, LNMF, LNMI, LRASCH 


REAL KAPPA 
REAL*4 ZMAT 
REAL*4 ZLEP 
REAL*4 SCHEME 


56 


no on no on n on on 


C 

DATA ZMAT /'MATS'/ 

DATA ZLEP /'LEAP'/ 

C 

DIMENSION TP ( 500) , TPSTAR( 500) , IMAXD2C 128) ,RH0C2) 
DIMENSION STOFFTC735) 

C 

LOGICAL FALSE 

INTEGER CORREC 

DATA Ml /-!/, FALSE/ *F. / 


COMMON/PRIME/ HPRIMEC 72 . 46,2) 
COMMON /PRIME/ UPRI MEC 7 2 , 46 , 2 ) 
COMMON/PRIME/ VPRIMEC 72 , 46 , 2 ) 


COMMON/UVH / U( 72, 46,2) 
COMMON/UVH / V(72,46, 2) 
COMMON/UVH/ HT(72, 46, 3) 
COMMON/UVH/ Q(72, 46, 2) 
COMMON/UVH/ HFORCEC 7 2,46,2) 
COMMON/UVH/ HH ( 72 , 46 ) 

COMMON/UVH/ UBAR (46,2) 
COMMON/UVH/ VBAR (46,2) 
COMMON/UVH/ HTBARC 2 ) 


COMMON/UVHP/ UP (2,2) 
COMMON/UVHP / VP ( 2 , 2 ) 
COMMON/UVHP/ HTP< 2,3) 


COMMON/DOT/ H(72,46,2,2) 

COMMON /DOT/ HU ( 72 , 46 , 2 , 2 ) 

COMMON /DOT/ HV ( 72 , 46 , 2 , 2 ) 

COMMON /DOT/ HQ (7 2, 46, 2) 

COMMON/DOT/ HD0TC72, 46) 
COMMON/DOT/ HUD0TC72, 46) 
COMMON/DOT/ HVDOT<72,46) 


COMMON /DOTP/ HP( 2, 2, 2) 

COMMON /DOTP/ HUP(2, 2, 2) 

COMMON /DOTP/ HVP(2, 2, 2) 

COMMON /DOTP/ HPD0TC2) 

COMMON /DOTP/ HUPDOTC 2 ) 
COMMON /DOTP/ HVPDOTC 2) 


COMMON /WORK/ SPACE ( 7 2 , 46 , 1 6 ) 
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C 


COMMON/PARM/ 

$ 

$ 

$ 

$ 

$ 

$ 

COMMON /DAMP/ 
COMMON /DAMP/ 


I MAX, JMAX,DL,DP,A,DT,K,KMAX,TITLEC20) , 
CON ( 46), 

COSPC 46), 

SI NP( 46), 

COSL( 72), 

SINLC 72), 

F< 46) 

DAMPG (46,36) 

DAMPA (46,36) 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


HT1XC72, 1), 
HT1YC72, 1), 
HU1XC72, 1), 
HV1YC72, 1), 
HUU1XC72, 1), 
HUV1XC72, 1), 
HVU 1 Y ( 72 , 1), 
HVV1YC72, 1), 


HT2XC72, 1 ) 
HT2YC72, 1 ) 
HU2XC72, 1 ) 
HV2YC72, 1 ) 
HUU2X (72, 1 ) 
HUV2XC 72 , 1 ) 
HVU2 Y (72, 1 ) 
HVV2YC72, 1 ) 


EQUIVALENCE 

c 

HT1XC 

EQUIVALENCE 

c 

HT2XC 

EQUIVALENCE 

c 

HT1YC 

EQUIVALENCE 

c 

HT2YC 

EQUIVALENCE 

c 

HU1XC 

EQUIVALENCE 

c 

HU2XC 

EQUIVALENCE 

c 

HV1 YC 

EQUIVALENCE 

c 

HV2YC 

EQUIVALENCE 

c 

HUU1XC 

EQUIVALENCE 

c 

HUU2XC 

EQUIVALENCE 

c 

HUV1XC 

EQUIVALENCE 

c 

HUV2XC 

EQUIVALENCE 

c 

HVU 1 Y C 

EQUIVALENCE 

c 

HVU2YC 

EQUIVALENCE 

c 

HVV1 YC 

EQUIVALENCE 

c 

HVV2YC 

IMAX =72 



JMAX = 46 



NMAX = IMAX/2 



1. 1 ) , SPACEC 1,1,01) ) 

1 . 1 ) , SPACEC 1,1,02) ) 

1. 1) , SPACEC 1, 1,03) ) 

1.1) .SPACEC 1 . 1 .04) ) 
1, 1),SPACEC1, 1,05) ) 

1.1) , SPACEC 1,1,06) ) 
1, 1 ) , SPACE Cl, 1,07) ) 

1. 1) , SPACEC 1, 1,08) ) 

1. 1 ) , SPACEC 1 , 1,09) ) 

1 , 1 ), SPACEC 1,1,10) ) 

1 , 1 ), SPACEC 1,1,11) ) 
1, 1), SPACEC 1,1,12) ) 

1. 1) , SPACEC 1, 1, 13) ) 
1, 1), SPACEC 1,1,14) ) 

1. 1) , SPACEC 1, 1, 15) ) 

1. 1) , SPACEC 1, 1, 16) ) 


CALL CLOCK 
C 

READC5, TSCHEM) 
READC5, INPUT) 
READC5, LEVELS) 
READC5,FILT) 
READC5, STRESS) 
READC 5 , FORCE) 
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c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


********************************************************* 
**** PHYSICAL PARAMETERS **** 

********************************************************* 


PI = 4 . 0*ATAN < 1 . OEO) 

G =9.81 

OMEGA = 0.7272205E-4 

A = 6.372E06 


********************************************************* 

**** TIME PARAMETERS **** 
**** DT IS IN SECONDS **** 
**** TMAX AND TPRT ARE IN HOURS **** 
********************************************************* 


TSHAP = TSHAP*3600 . 


I F< SCHEME. EQ.ZMAT) MATS = 1 
I F< SCHEME . EQ . ZLEP) MATS = 0 


CHANGE TIME UNITS TO SECS. 

TM = TMAX 

TO = TSTART 

TMAX = TMAX * 3600.0 
DO 40 1=1,500 
TPSTAR(I) = I *TPRT 

40 TPC I ) = ( TPSTARt I ) +T0) * 3600.0 


C 

C ********************************************************* 
C **** DEFINE DOMAIN LIMITS **** 

C ********************************************************* 
C 

JSP = 1 
JNP = JMAX 


DP = PI / ( JNP-JSP ) 
DL = 2 . 0*PI /IMAX 
DPSTAR = DP* 180/PI 


C 

C ********************************************************* 
C **** BEGIN MAIN INTEGRATION LOOP **** 

C ********************************************************* 
C 

15 CONTINUE 
C 

C ********************************************************* 
C **** CHOOSE PREDICTOR OR CORRECTOR FOR MATSUNO **** 
C ********************************************************* 
C 

I F< MATS . EQ . 1 ) CORREC = CORREC * Ml 
C 

DO 2 K= 1 , KMAX 
C 

C ********************************************************* 
C **** SET CURRENT VALUES **** 

C ********************************************************* 
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DO 101 J=JSP, JNP 
DO 101 1=1, IMAX 

HH< I , J) = RH0<K>*< HTd , J, 1)-HT<I,J,2> ) + HT(I,J,2) 

H<I,J,K,N) = HTCI,J,K) - HT( I , J,K+ 1 ) 

HUd,J,K,N> = Hd , J,K,N) * U< I , J,K) 

HV<I,J,K,N) = Hd , J,K,N) * Vd , J,K) 

101 CONTINUE 

DO 102 M= 1 , 2 

HP(M,K,N) = HTP(M,K) - HTP(M,K+1) 

HUPCM,K,N) = HP(M,K,N) * UP(M,K> 

HVPCM,K,N) = HP(M,K,N) * VP(M,K> 

102 CONTINUE 

********************************************************* 

**** CALL TRACER ROUTINE **** 

********************************************************* 

IF(LTRACE) CALL TRACER ( H(1,1,K,N), U(1,1,K), V(1,1,K), 

$ Q( 1 , 1 ,K) , HQd,l,K), IMAX, JMAX ) 

********************************************************* 

**** CALCULATE AVERAGE QUANTITIES **** 

********************************************************* 

DO 110 J=JSP, JNPM2 
JP1=J+1 
JP2=J+2 

I = IMAX- 1 
I P 1 = IMAX 

DO 120 IP2= 1 , IMAX 

U1X 
V1X 

HUIXd , J) 

HUUIXd, J) 

HUV1XCI , J) 

HTlXd , J) 

C 

U2X 
V2X 

HU2X< IP1 , J) 

HUU2XC IP1 , J) 

HUV2X< IP 1 , J) 

HT2X( I P 1 , J) 

C 

U 1 Y 
V 1 Y 

HV 1 Yd , J) 

HVUlYd, J) 

HVV1YCI, J) 

HT 1 Y< I , J) 


= UdPl , J,K) + U(I,J,K) 

= V< IP1 , J, K> + Vd , J , K ) 

= HUdP 1 , J, K, N> + HUCI , J,K,N) 
= HU 1XC I , J) * U1X 
= HU1X(I,J) * V1X 
= HHdPl , J) + HHd , J) 

= U ( IP2 , J, K) + Ud , J,K) 

= V C I P2 , J, K ) + Vd,J,K) 

= HUCIP2, J ,K,N) + HUd , J, K,N) 
= HU2X( IP1 , J) * U2X 
= HU2X(IP1,J) * V2X 
= HHCIP2, J) + HHd , J) 


= U(I , JP1 ,K) + Ud , J, K) 

= Vd , JP1 , K> + Vd , J, K) 

= HVd , JP1 ,K,N) *COSP( JP1) + HVC I , J,K,N>*COSP< J) 
= HV1YCI , J) * U1Y 
= HV 1 Yd , J) * V 1 Y 
= HHd , JP1 ) + HH d , J ) 
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c 

U2Y = UCI,JP2,K> + U(I,J,K> 

V2Y = VCI,JP2,K) + VCI,J,K> 

HV2 YC I , JP 1 ) = HV< I , JP2,K,N)*C0SPC JP2) + HV( I , J,K,N)*COSPC J) 
HVU2 YC I , JP 1 ) = HV2YC I , JP 1 ) * U2Y 
HVV2Y C I , JP 1 ) = HV2YC I , JP 1 ) * V2Y 
HT2 Y( I , JP 1 ) = HHC I , JP2) + HHCI,J> 

C 

I = IP1 
I P 1 = I P2 
120 CONTINUE 
110 CONTINUE 

********************************************************* 

**** CALCULATE AVERAGE QUANTITIES NEAR POLE **** 

********************************************************* 

J = JNPM 1 
JP1=J+1 
C 

I =IMAX- 1 
I P 1 = I MAX 

DO 130 IP2= 1 , IMAX 
C 

U1X = UCIP1,J,K> + U(I,J,K> 

V1X = V(IP1,J,K) + VCI,J,K) 

HU IX ( I , J) = HUCIP1,J,K,N) + HU( I , J,K,N) 

HUU1XC I , J) = HU1X< I , J) * U1X 

HUV 1 X ( I , J ) = HU1XCI,J> * V1X 

HT1X< I , J) = HH< IP 1 , J) + HH ( I , J ) 

C 

U2X = U(IP2,J,K) + UC I , J, K) 

V2X = V(IP2,J,K) + V(I,J,K) 

HU2X< IP1 , J) = HU(IP2,J,K,N) + HU(I,J,K,N) 

HUU2X( IP 1 , J) = HU2X< IP1 , J) * U2X 
HUV2X( IP 1 , J ) = HU2X( IP1 , J) * V2X 
HT2X( IP1 , J) = HHC IP2, J) + HH(I,J) 

C 

U 1 Y = UC I , JP 1 , K ) + UCI,J,K) 

V 1 Y = VC I , JP 1 , K ) + VCI,J,K> 

HV 1 Y C I , J ) = HV CI,JP1,K,N) *COSP C JP 1 ) + HV C I , J , K , N ) *COSP C J ) 

HVU 1 Y C I , J ) = HV 1 YC I , J) * U1Y 

HVV 1 Y C I , J ) = HV 1 YC I , J) * V 1 Y 

HT1YCI.J) = HHC I , JP1 ) + HHC I , J) 

C 

HV2YC I , JSP) = 0.0 

HVU2YC I , JSP) = 0.0 

HVV2YC I , JSP) = 0.0 

HT2YC I , JSP) = HHCI+IMAXD2C I ) , JSPP1 ) + HHCI,JSPP1) 

C 

HV2YC I , JNP) = 0.0 

HVU2YC I , JNP) = 0.0 

HVV2YC I , JNP) = 0.0 

HT2Y C I , JNP) = HHC I +IMAXD2C I ) , JNPM 1 ) + HHCI,JNPM1> 
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c 


I = I p 1 
IP 1 = I P2 
130 CONTINUE 
C 

C ********************************************************* 

c **** CALCULATE TIME TENDENCIES **** 

C ********************************************************* 

c 

DO 150 J=JSPP1 , JNPM1 
JP1=J+1 
JM1 =J- 1 
C 


I = IMAX 

IM 1 = IMAX- 1 

DO 160 I P 1 = 1 , I MAX 


FSTAR = F(J) + U(I,J,K)*SINP(J) * CON(J) 

PHIZ = G * H<I,J,K,N) 


C 

C 


c 

c 


c 

c 


c 


HUCOR = 

FSTAR 

* HVCI,J,K,N) 


HVCOR = 

FSTAR 

* HU( I , J, K, N) 


DHTDL = 

( HT1X( I , J) -HT 1X< IM 1 , J) ) *DL24 


$ 

“ 

( HT2XCIP1,J)-HT2X(IM1,J) 

) *DL4 

DHTDP = 

C HT 1 Y ( I , J) -HT 1 Y < I , JM 1 ) ) *DP24 


$ 


< HT2YCI ,JP1 )-HT2Y( I , JM 1 ) 

) *DP4 

HUPRES = 

PHIZ 

* CON(J) * DHTDL 


HVPRES = 

PHIZ 

* AINV * DHTDP 



HDOTCI, J) 


= -CON(J) * ( 



$ 


( 

HU 1 X ( 

I , J ) - HU 1 X ( I M 1 , J ) 

) *DL24 

$ 


- < 

HU2X< IP1 , J ) - HU2XC IM 1 . J ) 

>*DL4 

$ 

+ 

( 

HV 1 Y( 

I , J ) - HV 1 Y ( I , JM 1 ) 

)*DP24 

$ 


- ( 

HV2 Y ( 

I ,JP1) - HV2 Y ( I ,JM1> 

) *DP4 ) 

HUDOT(I, J) 


= -CON(J) * ( 



$ 


( 

HUU 1 X( 

I , J ) - HUU 1 X ( I M 1 , J ) 

) *DL44 

$ 


" ( 

HUU2X( IP 1 , J ) - HUU2X ( IM 1 , J ) 

) *DL8 

$ 

+ 

( 

HVU 1 Y < 

I , J ) - HVU 1 Y C I , JM 1 ) 

) *DP4 4 

$ 


“ < 

HVU2YC 

I ,JP1) - HVU2 Y ( I , JM 1 ) 

) *DP8 ) 

$ 

+ 

HUCOR 

- HUPRES 



HVDOT ( I , J) 

= -CON(J) * 

< 




$ 

( 

HUV 1X( 

I 

, J 

) - HUV1XCIM1, J ) 

)*DL44 

$ 

“ ( 

HUV2X( IP 1 

, J 

) - HUV2XCIM1, J ) 

) *DL8 

$ 

+ ( 

HVV1YC 

I 

, J 

) - HVV 1 Y ( I , JM 1 ) 

) *DP44 

$ 

- ( 

HVV2Y< 

I 

,JP1) - HVV2Y ( I , JM 1 ) 

) *DP8 ) 

$ 

- HVCOR 

- HVPRES 




I M 1 = I 
I =1 P 1 
160 CONTINUE 
150 CONTINUE 
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********************************************************* 
**** CALCULATE TIME TENDENCIES AT POLES **** 

********************************************************* 


C 

C 

C 


c 


c 

c 


DO 170 M= 1 , 2 
MSGN = < - 1 ) **M 

JSTAR1 = JSP+1 + (M-l )*( JNP-1 - (JSP+l) ) 

JSTAR2 = JSP+2 + (M-l )*( JNP-2 - (JSP+2) ) 

JPOLE = JSP + CM— 1 > * < JNP - JSP ) 

SUM 11 = 0.0 
SUM 12 = 0.0 
SUM21 = 0.0 
SUM22 = 0.0 
SUM41 = 0.0 
SUM42 = 0.0 
SUM31 = 0.0 
SUM32 = 0.0 
SUM51 = 0.0 
SUM52 = 0.0 
DO 180 1=1, IMAX 

SUM 1 1 = SUM 1 1 + HV< I , JSTAR 1 , K , N ) 

SUM 1 2 = SUM 1 2 + HV( I , JSTAR2 , K , N ) 

SUM21 = SUM21 + HV< I , JSTAR 1 , K, N) * ( -U ( I , JSTAR 1 , K) *SI NL( I ) 

$ - MSGN * V( I , JSTAR 1 ,K)*COSL( I ) ) 

SUM22 = SUM22 + HV ( I , JSTAR2 , K , N ) * ( -U ( I , JSTAR2 , K ) *SINL( I ) 

$ - MSGN * V(I ,JSTAR2,K)*C0SL(I ) ) 

SUM41 = SUM41 + HV ( I , JSTAR 1,K,N)*( MSGN * U ( I , JSTAR 1 , K ) *COSL( I ) 
$ - V( I , JSTAR1 ,K)*SINL( I ) ) 

SUM42 = SUM42 + HV( I , JSTAR 2 , K , N) *( MSGN * U( I , JSTAR2 , K) *COSL( I ) 
$ - V(I,JSTAR2,K)*SINL(I) ) 

SUM31 = SUM31 + COSL(I) * HH(I,JSTAR1) 

SUM32 = SUM32 + COSL(I) * HH(I,JSTAR2) 

SUM51 = SUM51 + SINL(I) * HH(I,JSTAR1> 

SUM52 = SUM52 + SINL(I) * HH(I,JSTAR2) 

180 CONTINUE 


C 


HPDOT(M) = MSGN * ( C11*SUM11 - C12*SUM12 ) 


C 

HUPDOT(M) = ( MSGN * ( C11*SUM21 - C12*SUM22 ) 

$ - G* HP(M,K,N)* ( C 1 1 *SUM3 1 - C12*SUM32 ) ) 

$ + F ( JPOLE) * HVP(M,K,N) 

C 

HVPDOT(M) = MSGN * ( ( C11*SUM41 - C12*SUM42 ) 

$ - G* HP( M, K, N ) * < C 1 1 *SUM5 1 - C12*SUM52 ) ) 

$ - F( JPOLE) * HUP(M,K,N) 

C 

170 CONTINUE 
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********************************************************* 
**** APPLY FOURIER FILTER TO TIME TENDENCIES **** 

********************************************************* 

IF< .NOT.HFILT) GOTO 325 

CALL FILTER < HDOT, IMAX, JMAX,STOFFT) 

CALL FILTER (HUDOT, IMAX, JMAX,STOFFT) 

CALL FILTER (HVDOT, IMAX, JMAX, STOFFT) 

325 CONTINUE 


********************************************************* 
**** ADD FORCING TERMS **** 

********************************************************* 

IF( .NOT.LFORCE) GOTO 260 

DO 265 J=JSPP1 , JNPM1 
DO 265 1=1, IMAX 

HDOT( I , J) = HDOT( I , J) + ALPHA*HFORCE< I , J, K) 

265 CONTINUE 

260 CONTINUE 

********************************************************* 
**** SUBTRACT INITIAL TENDENCY FOR PERTURBATION EXP. **** 
********************************************************* 

IFCMATS.NE. 1 ) GOTO 285 
IF( . NOT . LPRIME) GOTO 285 
I F( CORREC . EQ . - 1 ) GOTO 285 

DO 275 J=JSPP1 , JNPM1 
DO 275 1=1, IMAX 

HDOT ( I , J ) = HDOT ( I , J ) - HPRIME( I , J, K) 

HUDOT ( I , J ) = HUDOT ( I , J ) - UPRIMEC I , J, K) 

HVDOT ( I , J) = HVDOT ( I , J) - VPRIMEC I , J, K) 

275 CONTINUE 

285 CONTINUE 

********************************************************* 
**** USE LEAP FROG OR MATSUNO TO ADVANCE SOLUTION **** 
********************************************************* 

Z = 2.0 

I FCNCALL . EQ . 1) Z = 1 .0 

ZDT = Z*DT 


64 



ooono o o ooooo ooooo 


c 


c 


DO 270 J=JSPP1 , JNPM1 
DO 270 1 = 1 , 1 MAX 
H(I,J,K,NP1) = H(I,J,K,NM1) 
HU<I , J,K,NP1 ) = HU(I,J,K,NM1) 
HV<I,J,K,NP1) = HV<I,J,K,NM1) 
270 CONTINUE 


+ ZDT* HDOTC I , J) 
+ ZDT*HUD0T< I , J) 
+ ZDT*HVD0T(I,J> 


DO 280 M= 1 , 2 
HP(M,K,NP1 ) 
HUP(M,K,NP1 ) 
HVP(M,K,NP1) 
280 CONTINUE 


HP<M,K,NM1) 

HUP(M,K,NM1) 

HVP(M,K,NM1) 


+ ZDT* HPDOT(M) 
+ ZDT*HUPDOT(M) 
+ ZDT*HVPDOT(M) 


********************************************************* 
**** END K-LOOP **** 

********************************************************* 

2 CONTINUE 

********************************************************* 
**** CALCULATE NEW VELOCITIES AND HEIGHTS **** 

********************************************************* 

DO 295 KK= 1 , KMAX 
K = KMAX - KK + 1 


DO 290 J=JSPP1 , JNPM1 
DO 290 1 = 1 ,1 MAX 
U(I,J,K) = 0.0 
V< I , J , K ) = 0.0 

HT( I , J , K ) = HT( I , J,K+ 1 ) + H(I,J,K, NP 1 ) 

IF(H<I,J,K,NP1) .LE.O.O) GOTO 290 
U(I,J,K) = HU < I , J , K, NP 1 ) / H(I ,J,K, NP 1 ) 

V( I , J , K) = HV(I , J,K,NP1 ) / H(I,J,K, NP 1 ) 

QCI , J,K) = HQ( I , J,K) / H(I ,J,K, NP 1 ) 

290 CONTINUE 

DO 300 M= 1 , 2 
UP(M,K) = 0.0 
VP(M,K) = 0.0 

HTP(M, K) = HTP(M, K+ 1 ) + HP(M,K,NP1) 

IF(HP(M,K,NP1 ) .LE.O.O) GOTO 300 
UPCM.K) = HUP(M,K,NP1 ) / HP(M.K.NPl) 

VP ( M , K ) = HVP ( M , K, NP 1 ) / HP(M,K,NP1) 

300 CONTINUE 

********************************************************* 
**** APPLY BOUNDARY CONDITIONS AT POLES **** 

********************************************************* 

DO 310 M= 1 , 2 
C 

MSGN = < - 1 ) **M 

JPOLE = JSP + < M — 1 > * < JNP - JSP ) 
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DO 320 1=1, IMAX 

U ( I , JPOLE, K ) = -UP(M,K)*SINL<I) + MSGN * VP(M,K) *C0SL( I ) 
V( I , JP0LE,K) = -MSGN * UP( M , K ) *C0SL< I ) - VP(M,K)*SINL( I ) 
HT(I, JPOLE, K) = HTP(M,K) 

Q(I, JPOLE, K) = HQ ( I , JPOLE, K)/HTP<M,K) 

320 CONTINUE 
C 

310 CONTINUE 
295 CONTINUE 

******************************* ************************** 

**** DO BACKWARD STEP FOR MATSUNO **** 

********************************************************* 

I F ( NCALL . NE . 1 ) GOTO 412 

N = NP1 

NSTEP = NSTEP + 1 
IF< NSTEP . EQ . 1 ) GOTO 15 
NSTEP = 0 

412 CONTINUE 

********************************************************* 

**** ADD FRICTION AND ZONAL FORCING **** 

********************************************************* 

IF( .NOT.LFRICT) GOTO 221 
DO 224 K= 1 , KMAX 

CALL FRICT ( U ( 1 , 1 , K) , V ( 1 , 1 , K ) , UBAR ( 1 , K ) , VBAR < 1 , K ) , 

$ IMAX, JMAX,DT,TAO) 

224 CONTINUE 

221 CONTINUE 

********************************************************* 
**** UPDATE TIME **** 

********************************************************* 

T = T+DT 

******** C************************************************* 

**** APPLY NORMAL MODE FILTER **** 

********************************************************* 

IF< . NOT . LNMF) GOTO 321 

DO 324 K= 1 , KMAX 

CALL NORMAL ( IK 1 , 1 ,K> ,V< 1 , 1 ,K) ,HT( 1 , 1 ,K> , 

$ HTBAR(K), IMAX, JMAX, DT, LNMI , LNMF, LRASCH > 

324 CONTINUE 
C 

321 CONTINUE 
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c ********************************************************* 

c **** APPLY SHAPIRO FILTER **** 

C ********************************************************* 

C 

IF( .NOT.GFILT) .GOTO 330 

IF< (T-TSTART) . LT . NEXT*TSHAP ) GOTO 330 

C 

DO 600 K= 1 , KMAX 

CALL SHAP ( U(1,1,K), IMAX, JMAX, 1) 

CALL SHAP ( V( 1 , 1 ,K) , IMAX, JMAX, 1 ) 

CALL SHAP (HT( 1 , 1 ,K) , IMAX, JMAX, 2) 

DO 610 M= 1 , 2 

JPOLE = JSP + (M- 1 )*( JNP-JSP) 

SI = 0. 

DO 620 1=1, IMAX 

SI = SI + HT(I, JPOLE, K) 

620 CONTINUE 

HTSTAR = SI /IMAX 
DO 630 1=1, IMAX 
HT<I, JPOLE, K) = HTSTAR 
630 CONTINUE 
610 CONTINUE 
600 CONTINUE 
C 

TSTAR = T/3600. 

WRITE (6,3322) TSTAR, NCALL , NTOT 
3322 FORMAT( IX, ’TIME : ’ , IX, F8 . 3, 1 X, ’ HRS . ' , 6X, ’ NCALL :',I3,4X, 
$ ’NTOT : ’ , I 4, 4X, 'SHAPIRO FILTER APPLIED’) 

I F (SCHEME . EQ . ZMAT ) NEXT = NEXT + 1 

I F( SCHEME . EQ . ZMAT) GOTO 330 

ICOUNT = ICOUNT+1 
IF( I COUNT . EQ . 1 ) GOTO 330 

ICOUNT = 0 
NEXT = NEXT + 1 
330 CONTINUE 

********************************************************* 
**** TEST FOR INSTABILITY ON HEIGHT FIELD **** 

********************************************************* 

DO 7 K= 1 , KMAX 
C 

HBAR = 0. 1 *HTBAR(K) 

DO 450 J=JSP, JNP 
DO 450 1=1, IMAX 

IF(HT( I , J,K) .LE. HBAR) FLAG = 1.0 
450 CONTINUE 
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C 

IFCFLAG.EQ.O.O) GOTO 7 

TSTAR = T/3600.0 

WRI TEC 6 , 480) TSTAR 

WRITE< 18,41 17) TITLE, K, TSTAR 

CALL HARM (HT(1,1,K), IMAX, JMAX, 1 8 ) 

4117 F0RMATC/1X,20A4, /IX, 'LEVEL ’,12,' BECAME UNSTABLE AT T 
$ F7 . 3, ' HOURS’,//) 

GOTO 490 
C 

7 CONTINUE 


********************************************************* 
**** PERFORM NORMAL MODE INITIALIZATION **** 

********************************************************* 


IF(.NOT.LNMI) GOTO 457 
DO 458 K= 1 , KMAX 

CALL NORMAL C U( 1 , 1 ,K) , V< 1 , 1 ,K) ,HTC 1 , 1 ,K> , 

$ HTBARCK), IMAX, JMAX, DT, LNMI , LNMF, LRASCH ) 

458 CONTINUE 
CALL POLE 

T = T-DT 

NMICNT = NMICNT + 1 
IFCNMICNT.LT.NMIMAX) GOTO 14 
LNMI = FALSE 

WRI TE( 45) TITLE, IMAX, JMAX 

WRI TEC 45) < CHTC I , J, KTPG) , I = 1 , IMAX) , J= 1 , JMAX) 

DO 459 K= 1 , KMAX 
WRITEC 45) K, T 

WRI TEC 45 ) CC QC I , J, K) , I = 1 , IMAX) , J= 1 , JMAX ) 

WRI TEC 45 ) CC U C I , J , K) , I = 1 , 1 MAX ) , J= 1 , JMAX ) 

WRITEC45) CC V C I , J, K> , I = 1 , I MAX ) , J= 1 , JMAX > 

WRITEC 45) C CHTC I , J,K), 1 = 1 , IMAX) , J= 1 , JMAX) 

459 CONTINUE 
CALL CLOCK 
GOTO 14 

457 CONTINUE 


********************************************************* 
**** PRINT OUTPUT **** 

********************************************************* 


IFCT-TPCIPN) .LT. 0.0) GOTO 500 
C 

490 CONTINUE 

DO 8 K= 1 , KMAX 
TSTAR=T/ 3600 . 0 
C 


WRITEC 45) 
WRITEC 45 ) 
WRITEC 45 ) 
WRI TEC 45 ) 
WRITEC 45 ) 


K, TSTAR 

CC QCI,J,K), 1=1, IMAX), J=l, JMAX) 
CC UCI,J,K), 1=1, IMAX), J=l, JMAX) 
CC VCI,J,K), 1=1, IMAX), J=l, JMAX) 
CCHTCI,J,K),I=1,I MAX ) , J= 1 , JMAX ) 
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C 

WRITE (6,749) TSTAR, NCALL, NTOT, K 
C 

8 CONTINUE 


********************************************************* 
**** CHECK FOR END OF RUN **** 

********************************************************* 

IFCFLAG . EQ. 1.0) STOP 
I PN= I PN+ 1 

500 CONTINUE 

TDAY = 86400.0 * AI NT( T/ 86 400 . 0 ) 

IF(T.EQ.TDAY) CALL CLOCK 

IF(T.GE.TMAX) CALL CLOCK 

IF(T.GE.TMAX) STOP 

********************************************************* 
**** SWITCH NM 1 , N, NP1 ARRAYS **** 

**************************** *********** ****************** 

IF(SCHEME.EQ.ZLEP) NCALL = NCALL+1 
GOTOC 431,432,433), NCALL 


C 


C 


C 


COMMON/PARM/ 


$ 

$ 

$ 

$ 


IMAX, JMAX,DL,DP, A,DT,K,KMAX, TITLE (20) , 
C0N( 46), 

COSPC 46), 

SINP< 46), 

COSL( 72), 


$ 


SI NE< 38), 


KTP = KMAX+1 

READ< 45,END=100) TITLE, IMAX, JMAX 

READ ( 45 , END= 100) ( < HT C I , J , KTP ) , I = 1 , 1 MAX ) , J= 1 , JMAX ) 


DO 10 K= 1 , KMAX 
DO 20 J= 1 , JMAX 
DO 30 1=1, IMAX 
HPRIMEd ,J,K) = 0.0 
UPRIMEd , J,K) = 0.0 
VPRIME( I , J,K) = 0.0 
30 CONTINUE 
20 CONTINUE 
10 CONTINUE 


200 CONTINUE 

DO 40 L= 1 , KMAX 
READ( 45 , END= 100) 
READ C 45 , END= 100) 
READ( 4 5, END = 100) 
READ(45,END=100) 
40 CONTINUE 


K TSTAR 

(( Ud, J,K), 1 = 1, IMAX), J=l, JMAX) 
(( Vd, J,K), 1 = 1, IMAX), J=l, JMAX) 
<<HTd,J,K),I = l,I MAX ) , J = 1 , JMAX) 
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DO 50 KK= 1 , KMAX 
K = KTP-KK 
DO 60 J= 1 , JMAX 
DO 70 1=1, IMAX 

Hd,J,K,l> = HTd,J,K)-HTd,J,K+l) 

70 CONTINUE 
60 CONTINUE 
C 

DO 80 J= 1 , JMAX 
DO 90 1=1, IMAX 

HPRIME< I , J,K) = H<I,J,K, 1) - HPRIMEC I , J,K> 

UPRIMEd , J,K) = U(I,J,K)*H(I,J,K, 1) - UPRIMEd, J,K> 
VPRIME< I , J,K) = Vd,J,K)*Hd,J,K, 1> - VPRIMEC I , J,K) 
90 CONTINUE 
80 CONTINUE 
50 CONTINUE 
C 

GOTO 200 
100 CONTINUE 
C 

DO 110 K= 1 , KMAX 
DO 120 J= 1 , JMAX 
DO 130 1=1, IMAX 

HPRIMEd, J,K) = HPRIMEd,J,K)/DT 
UPRIMEd, J,K) = UPRIMEd, J,K)/DT 
VPRIME< I , J,K) = VPRIMEC I,J,K)/DT 
130 CONTINUE 
120 CONTINUE 
1 10 CONTINUE 
C 

REWIND 45 

RETURN 

END 

SUBROUTINE UVHR (TSTART,LFORCE,LNMI ) 

DIMENSION GEOP( 46,72), TOPOG< 72 , 46 ) 

REAL*8 XLAB(IO) 

433 NM1 = N 
N = NP 1 
NP1 = NM1 
C 

IF(NCALL-1 .LE.NLEAP) GOTO 439 
NSA = N 
N = NP1 
NP 1 = NSA 
NM 1 = N 
NCALL = 1 
GOTO 439 
C 

431 NSA = NM 1 
NM 1 = NP 1 
NP 1 = NSA 
N = NM 1 
GOTO 439 

C 

432 NP 1 = NM 1 
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c 

439 CONTINUE 
C 

NTOT =NTOT +1 
GO TO 15 


749 FORMAT< IX, 'TIME : ' , IX, F8 . 3, IX , ' HRS . ' , 6X , ' NCALL :',I3,4X, 

$ 'NTOT :', 14, 4X, 'LEVEL :',I2) 

480 F0RMAT(//,5X, 'HEIGHT FIELD BECAME UNSTABLE AT T=',G12.4,' HRS.') 
105 FORMATC IX, 'NO. OF GRID POINTS IN LONGITUDE: ',13,/, IX, 

$ ' NO . OF GRID POINTS IN LATITUDE: ’,13,/, IX, 

$'GRID SIZE: ',F5.2,' DEG X',F5.2,' DEG’,/, IX, 

$'TIME STEP: ’,F6.1,’ SECS.’,/, IX, 

$ ’ PRI NTI NG EVERY ' ,F6 . 2, IX, 'HR(S) UP TO ' ,F7 . 2, IX, 'HR(S) . STARTI ’ , 
$’NG AT ' ,F6. 1 , IX, 'HOURS' ) 

108 FORMAT( IX, 'DENSITY RATIO: ’,F5.2> 

109 FORMAT ( IX, 'NO. OF LEVELS: ’,13) 

END 

SUBROUTINE PRIMES 

COMMON/PRIME/ HPRIMEC72, 46, 2) 

COMMON/PRIME/ UPRIME(72, 46, 2) 

COMMON/PRIME/ VPRIME(72, 46, 2) 


COMMON/UVH/ U<72,46,2) 
COMMON/UVH/ V(72,46,2) 
COMMON/UVH/ HT<72, 46, 3) 
COMMON/UVH/ 0(72,46,2) 
COMMON/UVH/ HF0RCE(72, 46,2) 
COMMON/UVH/ HH (72,46) 

COMMON/UVH/ UBAR (46,2) 
COMMON/UVH/ VBAR (46,2) 
COMMON/UVH/ HTBAR( 2 ) 


COMMON/DOT/ H(72,46,2,2) 

COMMON/DOT/ HU ( 72 , 46 , 2 , 2 ) 

COMMON/DOT/ HV ( 72 , 46 , 2 , 2 ) 

COMMON/DOT/ HQ(72,46,2) 

COMMON/DOT/ HD0T(72,46) 
COMMON/DOT/ HUD0T(72,46) 
COMMON/DOT/ HVD0T(72.46) 
EQUIVALENCE ( XLAB( 1 ) , TITLE( 1 ) ) 
LOGICAL LFORCE, LNMI 
NAMELIST /LABEL/ XLAB 

COMMON/UVH/ U(72,46,2) 

COMMON/UVH/ V(72,46,2) 

COMMON/UVH/ HT(72, 46, 3) 

COMMON/UVH/ 0(72,46,2) 

COMMON/UVH/ HF0RCE(72, 46,2) 
COMMON/UVH/ HH(72, 46) 
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COMMON/UVH/ UBARC 46 , 2 ) 

COMMON / UVH / VBAR (46,2) 

COMMON/UVH/ HTBAR( 2) 

C 

COMMON/PARM/ I MAX, JMAX, DL, DP , A , DT, K, KMAX, TI TLE( 20) , 


$ 

CON< 

46) 

$ 

COSP< 

46) 

$ 

SI NP< 

46) 

$ 

COSL< 

72) 

$ 

SI NL( 

72) 

$ 

F ( 

46) 

OMEGA = 
G = 9.8 

0. 7272205E- 

4 


C 

KTP = KMAX+1 
C 

C ********************************************************* 
C **** READ IN INITIAL U,V,H AND TOPOGRAPHY **** 

C ********************************************************* 
C 
C 

READ( 46 , END = 100) TITLE, I MAX, JMAX 

READ (46,END=100) ( ( HT< I , J , KTP ) , I = 1 , 1 MAX ) , J= 1 , JMAX ) 
READC5, LABEL) 

200 CONTINUE 

DO 500 L= 1 , KMAX 

READ ( 46 , END = 1 00) K, TSTAR 

READ< 46,END= 100) (( Q(I,J,K),I = 1, I MAX) , J= 1 , JMAX) 

READ (46,END=100) (( U(I,J,K),I = 1, I MAX) , J= 1 , JMAX) 

READ ( 46 , END= 100) (( V ( I , J , K) , I = 1 , 1 MAX ) , J= 1 , JMAX ) 

READC 46,END=100) ( (HT( I , J,K) , 1 = 1 , 1 MAX ) , J= 1 , JMAX ) 

500 CONTINUE 
C 

IFCTSTAR.GE. TSTART) GOTO 100 
GOTO 200 
100 CONTINUE 
C 

C ********************************************************* 
C **** READ IN FORCING TERMS **** 

C *****************^****)*C*************>»f****)lc**5«C************ 
C 

IFC .NOT.LFORCE) GOTO 400 
DO 300 K= 1 , KMAX 

READC 49) ( (HFORCEC I , J,K) , 1 = 1 , I MAX) , J= 1 , JMAX) 

300 CONTINUE 
400 CONTINUE 
C 

C **********^***5»C*****************^^*********************** 
C **** WRITE INITIAL U,V,H AND TOPOGRAPHY **** 

C ********************************************************* 

c 

CALL POLE 
IF(LNMI) RETURN 

TSTART = TSTAR*3600 . 0 


o o 


WRITEC 45) TITLE, IMAX, JMAX 

WRITEC 45) CCHTCI,J,KTP),I=1, IMAX) , J= 1 , JMAX) 

DO 600 K= 1 , KMAX 

WRI TE( 45 ) K, TSTAR 

WRITEC45) (( Q(I,J,K), 1=1, IMAX), J= 1 , JMAX) 
WRITEC45) (( UCI,J,K), 1=1, IMAX), J= 1 , JMAX) 
WRI TEC 45 ) (C V(I,J,K), 1 = 1, IMAX), J=l, JMAX) 
WRI TEC 45) CCHTCI,J,K),I = 1,I MAX ) , J= 1 , JMAX ) 
600 CONTINUE 
C 

RETURN 

END 

SUBROUTINE FILTER CQ, IM, JM, STO) 

DIMENSION Q C I M , JM ) , QTEMPC 360) , ST0C735) 
COMMON /DAMP/ DAMPGC46,36) 

COMMON /DAMP/ DAMPAC46,36) 


NMAX = IM/2 
JMM 1 = JM- 1 
ZIM = l./FLOATCIM) 

C 

DO 10 J= 2 , JMM 1 
C 

DMIN = 1.0 
DO 15 N= 1 , NMAX 

DMIN = AMIN1 CDAMPGC J,N) ,DMIN) 

15 CONTINUE 

IFCDMIN.GT.0.99) GOTO 10 
C 

DO 20 1=1, IM 
QTEMPCI) = Q C I , J ) 

20 CONTINUE 
C 

CALL RFFTF C I M , QTEMP , STO ) 

DO 30 NWAVE= 1 , NMAX 
N2 = NWAVE*2 

QTEMPCN2) = QTEMPCN2)* DAMPGC J, NWAVE) 
IFCNWAVE .EQ. NMAX) GO TO 30 

QTEMPCN2 + 1 ) = QTEMP CN 2+ 1 ) *DAMPG C J , NWAVE) 
30 CONTINUE 

CALL RFFTB C I M , QTEMP, STO ) 

C 

DO 40 1=1, IM 
Q C I , J ) = QTEMPC I ) *ZIM 
40 CONTINUE 
C 

10 CONTINUE 
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c 


$ 

$ 

$ 

$ 

$ 

$ 


RETURN 

END 

SUBROUTINE COEFF (HTBAR) 

DIMENSION HTBAR( 2) 

COMMON /PARM / I MAX , JMAX , DL , DP , A , DT , K , KMAX , TI TLE( 20 ) 

CON ( 46), 

COSP< 46), 

SINP( 46), 

COSL( 72), 

SI NL( 72), 

F( 46) 

COMMON /DAMP/ DAMPG(46,36) 

COMMON /DAMP/ DAMPA<46,36) 


********************************************************* 
**** SET DEFAULT VALUES FOR DAMPING COEFFICIENTS **** 
********************************************************* 


C 


c 


c 


c 


NMAX = I MAX / 2 

GH = 9.81 *HTBAR< 1 ) 

CG = SQRT(GH) 

UBAR = 50.0 

DNORM = (DL*A) / (DT*(UBAR+CG) ) 

DO 10 J = 1 , JMAX 
DM IN = 1.0 
DO 20 N = 1 , NMAX 
THETA = N*DL 
THETA2 = THETA + THETA 

W = DNORM * 6.0 / ( 8 . *SI N ( THETA ) 
D = COSP(J)*W 
DAMPG(J,N) = AMIN 1(1. 0,D) 

DMIN = AMI N 1 <DAMPG( J, N ) , DMI N ) 
20 CONTINUE 

IF (DMIN .LT. 0.99) GOTO 10 
DO 30 N = 1 , NMAX 
DAMPG( J, N) = 1.0 
30 CONTINUE 
10 CONTINUE 

RETURN 

END 

SUBROUTINE POLE 

COMMON/UVH/ U<72,46,2) 

COMMON/UVH/ V(72,46,2) 

COMMON/UVH/ HT(72, 46, 3) 

COMMON/UVH/ Q<72, 46,2) 

COMMON/UVH/ HF0RCE(72, 46,2) 
COMMON/UVH/ HH(72, 46) 


- SIN(THETA2) ) 


COMMON/UVH/ UBAR( 46 , 2 ) 
COMMON/UVH/ VBARC 46,2) 
COMMON/UVH/ HTBAR( 2) 
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COMMON/UVHP/ UPC2, 2) 

COMMON/UVHP/ VP< 2, 2) 

COMMON/UVHP/ HTPC2,3) 

C 

COMMON /PARM / I MAX , JMAX , DL , DP , A , DT , K , KMAX , TI TLEC 20) , 

$ CONC 46), 

$ COSPC 46), 

$ SI NP< 46), 

$ COSL( 72), 

$ SINLC 72), 

$ FC 46) 

C 

OMEGA = 0. 7272205E-4 
G = 9.8 
C 

DO 50 M= 1 , 2 

JPOLE = 1 + CM- 1 )*C JMAX- 1 ) 

HTP ( M , KMAX+ 1 ) = HTC 1 , JPOLE, KMAX+1 ) 

DO 70 1=1,1 MAX 

HT(I, JPOLE, KMAX+1) = HTP(M, KMAX+ 1 ) 

70 CONTINUE 
50 CONTINUE 
C 

DO 80 K= 1 , KMAX 
DO 90 M= 1 , 2 

JPOLE = 1 + CM- 1 )*C JMAX- 1 ) 

MSGN = ( - 1 ) **M 
C 

UPOLE = 0.0 
VPOLE = 0.0 
DO 100 I =1,1 MAX 
SL = SINCCI-l) *DL) 

CL = COSC ( I - 1 ) *DL) 

UPOLE = UPOLE - U(I , JPOLE, K)*SL - MSGN * VC I , JPOLE, K)*CL 
VPOLE = VPOLE + MSGN * UC I , JPOLE, K) *CL - VC I , JPOLE, K)*SL 

100 CONTINUE 
C 

HTPCM,K) = HTC1, JPOLE, K) 

UPCM,K) = UPOLE/ IMAX 
VPCM,K) = VPOLE/ IMAX 
C 

DO 101 1=1, IMAX 
SL = SINCCI-1) *DL) 

CL = COSC CI-1 )*DL) 

UCI, JPOLE, K) = -UPCM,K)*SL + MSGN * VPCM,K)*CL 

VCI, JPOLE, K) = -MSGN * UPCM,K)*CL - VPCM,K)*SL 

HTCI, JPOLE, K) = HTPCM,K) 

101 CONTINUE 
C 

90 CONTINUE 
C 

80 CONTINUE 
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RETURN 

END 

SUBROUTINE AVERAG <Q, IM, JM, QBAR, LNMI ,LNMF) 

DIMENSION Q( IM, JM) , COSPC180) 

LOGICAL LNMI, LNMF 
DATA WR /9358 . 2996/ 

C 

PI = 4 . *ATAN( 1 . ) 

DL = 2 . *PI / IM 
DP = PI/CJM-l) 

DO 20 J=1,JM 

PHI = C J- 1 ) *DP - PI/2. 

COSP(J) = COS ( PH I ) 

20 CONTINUE 
C 

QSMJ = 0.0 
DO 1 10 J= 1 , JM 
QSMI = 0.0 
DO 120 1=1, IM 

QSMI = QSMI + Q<I,J> 

120 CONTINUE 

QSMJ = QSMJ + QSMI *COSP( J) 

1 10 CONTINUE 

QBAR = QSMJ * DL*DP/(4.*PI) 

C 

IF(. NOT. LNMI .AND. .NOT. LNMF) RETURN 
DO 130 J= 1 , JM 
DO 140 1=1, IM 

Q(I,J) = Q < I , J ) - QBAR + WR 
140 CONTINUE 
130 CONTINUE 
C 

QSMJ = 0.0 
DO 150 J= 1 , JM 
QSMI = 0.0 
DO 160 1=1, IM 

QSMI = QSMI + Q ( I , J) 

160 CONTINUE 

QSMJ = QSMJ + QSMI *COSP( J ) 

DLSTAR = DL* 180/PI 
KTPG = KMAX + 1 
C 

JNPM1 = JNP-1 
JNPM2 = JNP-2 
JSPP2 = JSP+2 
JSPP1 = JSP+1 
C 

C ********************************************************* 

C **** MODEL CONSTANTS **** 

C ********************************************************* 
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AINV = l./A 

DL2 = 1 . /<2.0*DL*3.0) 

DP2 = 1 . /(2.0*DP*3.0) 

DL4 = 1 . /(4.0*DL*3.0> 

DP4 = 1 . /<4.0*DP*3.0) 

DL8 = 1 . /(8.0*DL*3.0) 

DP8 = 1 . /<8.0*DP*3.0) 

DL24 = 4 . *DL2 
DP24 = 4 . *DP2 
DL44 = 4 . *DL4 
DP 4 4 = 4 . *DP4 
IDIM = IMAX 
JDIM = JMAX 
IMAXH = IMAX/2 

Cll = 4.* SIN< DP ) / ( 3 . *A*IMAX*( 1 • 0-C0S( DP >)> 

C 1 2 = SIN(2.*DP)/(3.*A*IMAX*( 1 . 0-C0S( 2 . *DP ) ) > 

CORREC = -1 
IF(MATS.EQ. 1 ) CORREC = 1 

RHO(l) = 1.0 
RHO< 2 ) = KAPPA 
DO 60 J=JSP, JNP 
PHI = ( J-JSP)*DP - PI/2.0 
COSP(J) = COS(PHI) 

SINP< J) = SIN(PHI) 

IF< J.NE. JSP .AND. J.NE.JNP) CON(J) = 1.0/ ( A*COSP< J) ) 
F ( J ) = 2 . 0*0MEGA*SINP( J) 

60 CONTINUE 

DO 70 1=1, IMAX 

IMAXD2( I ) = IMAXH - IMAX*< I /(IMAXH +1) ) 

COSL(I) = COS( < 1-1 )*DL) 

SINL(I) = SIN((I-1)*DL) 

70 CONTINUE 


********************************************************* 
**** READ INITIAL DATA **** 

********************************************************* 


IF(LPRIME) CALL PRIMES 

CALL UVHR ( TSTART, LFORCE, LNMI ) 

DO 55 K= 1 , KMAX 

CALL AVERAG (HT( 1 , 1 ,K) , IMAX, JMAX, HTBAR(K) , LNMI ,LNMF) 
CALL ZONAL ( U< 1 , 1 ,K> , IMAX, JMAX, UBAR< 1 ,K> ) 

55 CONTINUE 
C 

c ********************************************************* 

C **** WRITE HEADER **** 

C ********************************************************* 


C 


IF(SCHEME.EQ.ZLEP) 
I F< SCHEME. EQ.ZMAT) 
IF(GFILT) 

IF(HFILT) 

IF(LNMF) 


WRITE(6,50) TITLE 
WRITE (6, 54) 
WRITE(6,53) 
WRITE(6,5 1 ) 
WRITE<6,52) 
WRITEC6.56) 
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IF(LNMI ) 


WRI TEC 6 , 59 ) NMIMAX 
WRI TE< 6 ,109) KMAX 
WRI TE< 6 ,108) KAPPA 
WRI TEC 6 ,105) I MAX , JMAX, DLSTAR, 

DPSTAR, DT, TPRT, TM, TO 

50 FORMATC IX,//, IX, 'FOURTH ORDER 2-LEVEL SHALLOW WATER MODEL' 

$ ,/ , 1 X, 20A4, / ) 

51 FORMATC IX, 'GLOBAL SHAPIRO FILTER ') 

52 FORMATC IX, 'HIGH-LATITUDE FOURIER FILTER') 

53 FORMATC IX, 'MATSUNO TIME SCHEME') 

54 FORMATC IX, 'LEAPFROG TIME SCHEME’) 

56 FORMATC IX, 'HIGH-LATITUDE NORMAL MODE FILTER') 

77 FORMATC / IX, 'LEVEL ',12,' SCALE HEIGHT =',F13.4/) 

59 FORMATC IX, 'NON-LINEAR NORMAL MODE INITIALIZATION: ' 12 IX 

$ 'ITERATIONS') ’ ’ 

DO 76 K = 1 , KMAX 

WRI TEC 6,77) K, HTBARCK) 

76 CONTINUE 

********************************************************* 

**** PERFORM NORMAL MODE INITIALIZATION **** 

************************************ ********************* 


IFC . NOT . LNMI ) GOTO 57 
DO 58 K= 1 , KMAX 

CALL NORMAL C UC 1 , 1 ,K),VC 1 , 1 ,K) ,HTC 1 , 1 ,K), 

$ HTBARCK), IMAX, JMAX, DT, LNMI, LNMF,LRASCH ) 

58 CONTINUE 

57 CONTINUE 


************************************************** ******* 
**** TIME-SETP INITIALIZATION **** 

********************************************************* 


TAVE = 0.0 
T = TSTART 

NCALL = 1 
NTOT = 1 
IPN = 1 
NEXT = 1 


NSTEP = 0 
NMICNT = 0 
C 

CALL RFFTI C IMAX,STOFFT) 

C 

c ********************************************************* 

c **** CALCULATE DAMPING COEFFICIENTS **** 

C **** FOR FOURIER FILTERING **** 

C ********************************************************* 
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c 


CALL COEFF (HTBAR) 
CALL CLOCK 


********************************************************* 
**** SET INITIAL TIME INDEX **** 

********************************************************* 


14 CONTINUE 
NM1 = 1 
N = 1 
NP1 = 2 
150 CONTINUE 

QBAR = QSMJ * DL*DP/<4.*PI ) 

RETURN 

END 


R; T= 2 . 04 / 4 . 45 19:39:35 


79 


1. Report No. 2. Government Accession No. 

NASA TM-86227 

4. Title and Subtitle 

Documentation of the Goddard Laboratory for 
Atmospheres Fourth-Order Two-Layer Shallow Water 
Model 


7. Author(s) 

Lawrence L. Takacs, Compiler 


9. Performing Organization Name and Address 
Goddard Modelling and Simulation Branch 
Goddard Space Flight Center 
Greenbelt, Maryland 20771 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


BIBLIOGRAPHIC DATA SHEET 


3. Recipient's Catalog No. 


5. Report Date 

FEBRUARY 1986 


6. Performing Organization Code 
611 


8. Performing Organization Report No. 
85B0505 


10. Work Unit No. 


1 1 . Contract or Grant No. 


13. Type of Report and Period Covered 
Technical Memorandum 


14. Sponsoring Agency Code 


15. Supplementary Notes 


Lawrence L. Takacs: Sigma Data Services Corporation. 


16. Abstract 

This documentation describes the theory and numerical treatment used in the 2-level GLA 
fourth-order shallow water model. This model was designed to emulate the horizontal finite- 
differences used by the GLA Fourth-Order General Circulation Model (Kalnay et al., 1983) 
in addition to its grid structure, form of high-latitude and global filtering, and time- 
integration schemes. A user s guide is also provided instructing the user on how to create in- 
itial conditions, execute the model, and post-process the data history. 


17. Key Words (Selected by Author(s)) 

shallow water model, global filtering, 
2-layer, documentation, time-integration 
schemes, user’s guide, high-latitude filtering, 
horizontal finite-differences 


18. Distribution Statement 
unclassified - unlimited 


subject category 47 


19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of Pages I 22. Price 
unclassified unclassified 84 A05 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1986 









National Aeronautics and 
Space Administration 
Code NIT-4 

Washington, D.C. 
20546-0001 


BULK RATE 

POSTAGE & FEES PAID 

NASA Washington, DC 
Permit No. G-27 


Official Business 

Penalty for Private Use, $300 




POSTMASTER: 


If UndeUverable (Section 1 58 
Postal Manual) Do Not Return 



